{
 "metadata": {
  "name": "Labs 3-7"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### This document was written by Spencer Lyon"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Lab 3"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 1"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# Given in lab text, converted to funciton\n",
      "def create_sq_mat(k):\n",
      "    \"\"\"create k x k matrix of ascending ints\"\"\"\n",
      "    return [range(i, i+k) for i in range(0, k**2, k)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# list-based function\n",
      "def mat_mul(a, b):\n",
      "    \"\"\"\n",
      "    Multipliy matrices a and b\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    a, b : lists\n",
      "        The square matrices of the same size to be multiplied.\n",
      "        They should be represented as lists of lists.\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    res : lists\n",
      "        The result of the matrix multiplication\n",
      "\n",
      "    \"\"\"\n",
      "    rows_a = len(a)\n",
      "    cols_a = len(a[0])\n",
      "    rows_b = len(b)\n",
      "    cols_b = len(b[0])\n",
      "    \n",
      "    if cols_a != rows_b:\n",
      "        raise ValueError('Incompatiable dimensions')\n",
      "    \n",
      "    # Make placeholder list for result\n",
      "    res = [[0 for row in range(cols_b)] for col in range(rows_a)] \n",
      "    for i in range(rows_a):\n",
      "        for j in range(cols_b):\n",
      "            for k in range(rows_b):\n",
      "                res[i][j] += a[i][k] * b[k][j]\n",
      "    \n",
      "    return res"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = create_sq_mat(100)\n",
      "%timeit mat_mul(a, a)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 280 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "b = np.asarray(a)\n",
      "%timeit np.dot(b, b)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.13 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from time import time\n",
      "import pandas as pd\n",
      "\n",
      "k_vals = [5, 10, 50, 200, 400]\n",
      "time_np = []\n",
      "time_list = []\n",
      "\n",
      "for k in k_vals:\n",
      "    b = create_sq_mat(k)\n",
      "    c = np.arange(k ** 2).reshape(k, k)\n",
      "    \n",
      "    # Do numpy time\n",
      "    t1 = time()\n",
      "    np.dot(c, c)\n",
      "    time_np.append(time() - t1)\n",
      "    \n",
      "    # Do list time\n",
      "    t1 = time()\n",
      "    mat_mul(b, b)\n",
      "    time_list.append(time() - t1)\n",
      "\n",
      "times = pd.DataFrame([time_np, time_list], columns=k_vals, index=['NumPy', 'Lists']).T\n",
      "times.index.name = 'k'\n",
      "times"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>NumPy</th>\n",
        "      <th>Lists</th>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>k</th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>5  </th>\n",
        "      <td> 0.000022</td>\n",
        "      <td>  0.000108</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>10 </th>\n",
        "      <td> 0.000008</td>\n",
        "      <td>  0.000643</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>50 </th>\n",
        "      <td> 0.000145</td>\n",
        "      <td>  0.040745</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>200</th>\n",
        "      <td> 0.009628</td>\n",
        "      <td>  2.325591</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>400</th>\n",
        "      <td> 0.122684</td>\n",
        "      <td> 24.958129</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "        NumPy      Lists\n",
        "k                       \n",
        "5    0.000022   0.000108\n",
        "10   0.000008   0.000643\n",
        "50   0.000145   0.040745\n",
        "200  0.009628   2.325591\n",
        "400  0.122684  24.958129"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 2"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = np.random.randn(50, 1, 3, 1)\n",
      "b = np.random.randn(50, 1, 3, 1)\n",
      "c = np.random.randn(1, 42, 3, 7)\n",
      "d = np.random.randn(7)\n",
      "\n",
      "case_1 = a + b\n",
      "case_2 = a / c\n",
      "case_3 = c * d\n",
      "\n",
      "msg = 'in 1 shape: %s\\nin 2 shape: %s\\nout shape: %s'\n",
      "print('Case 1')\n",
      "print(msg % (a.shape, b.shape, case_1.shape))\n",
      "\n",
      "print('\\n\\nCase 2')\n",
      "print(msg % (a.shape, c.shape, case_2.shape))\n",
      "\n",
      "print('\\n\\nCase 3')\n",
      "print(msg % (c.shape, d.shape, case_3.shape))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Case 1\n",
        "in 1 shape: (50, 1, 3, 1)\n",
        "in 2 shape: (50, 1, 3, 1)\n",
        "out shape: (50, 1, 3, 1)\n",
        "\n",
        "\n",
        "Case 2\n",
        "in 1 shape: (50, 1, 3, 1)\n",
        "in 2 shape: (1, 42, 3, 7)\n",
        "out shape: (50, 42, 3, 7)\n",
        "\n",
        "\n",
        "Case 3\n",
        "in 1 shape: (1, 42, 3, 7)\n",
        "in 2 shape: (7,)\n",
        "out shape: (1, 42, 3, 7)\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Lab 4"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 1"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "x = np.linspace(-5, 5, 20)\n",
      "\n",
      "plt.plot(x, 3 * x, 'kD')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD9CAYAAACoXlzKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEwRJREFUeJzt3X9oVfUfx/HXbZtQaTnBLXOjq06bk5lXZlpgXchNTaeL\nytIwuCuRhKAsEv+opuCPfkiYYvmHjSQo8485iBwr4ppEtswJklSDNvw1TXNGZTSdn+8fsfvdDze3\n3XPuuedzng8Y7N7de8/7QL04fe7nvAoZY4wAANa4yesBAADOItgBwDIEOwBYhmAHAMsQ7ABgGYId\nACyTVLBXVlYqNzdXxcXFieeqqqqUl5enSCSiSCSiurq6pIcEAAxcUsEei8V6BXcoFNLq1avV2Nio\nxsZGzZs3L6kBAQCDk1Swz549W9nZ2b2e554nAPBOphsfum3bNu3evVslJSXasmWLRo4c2e3voVDI\njcMCgPUGcuHs+Jenzz33nJqbm3X06FGNGTNGL730Up/D2frz+uuvez4D58f5BfH83D63ixcvqqSk\n5LqZVlJSoosXL7p6/IFyPNhzcnIUCoUUCoX07LPPqqGhwelDAIAnsrOzVV9f3yvcS0pKVF9ff92l\naS84Huytra2J32tqarrtmAEAv+sZ7ukW6lKSa+xLly7VgQMHdOHCBeXn52vdunWKx+M6evSoQqGQ\nxo0bp507dzo1q29Eo1GvR3AV5+dvNp9fqs6tM9xjsZiqq6vTKtQlKWQGs3Dj1EFDoUGtFwEABp6d\n3HkKAJYh2AEEWltbmyoqKtTW1ub1KI4h2AEEVltbm8rKylRbW6uysjJrwp1gBxBInaF++PBhSdLh\nw4etCXeCHUDg9Az1TraEO8EOIHBisVivUO90+PBhxWKxFE/kLLY7Agicvq7YpfS84agT2x0BoA9+\nqQYYKoIdQCD5oRpgqFiKARBobW1taVsN0NNAs5NgBwCfYI0dAAKKYAcAyxDsAHzPxr6XZBDsAHzN\n1r6XZBDsAHzL5r6XZBDsAHzJ9r6XZBDsAHzJ9r6XZLCPHYAv+bXvJRnsYwdgNdv7XpJBsAPwLZv7\nXpLBUgwA3/NT30sy6IoBAMuwxg4AAUWwA4BlCHYAaYG+F+cQ7AA8R9+Lswh2AJ6i78V5SQV7ZWWl\ncnNzVVxcnHju4sWLKi0t1aRJk1RWVqZLly4lPSQAO9H34o6kgj0Wi6murq7bc5s3b1Zpaal++eUX\nPfTQQ9q8eXNSAwKwF30v7kh6H3tLS4vKy8t17NgxSVJhYaEOHDig3NxcnT17VtFoVD/99FP3g7KP\nHYCC2feSjIFmZ6bTBz537pxyc3MlSbm5uTp37tx1X1dVVZX4PRqNKhqNOj0KgDTXWQnQM9wJ9f/E\n43HF4/FBv8/xK/bs7Oxu62KjRo3SxYsXux+UK3YAXXS9cifU++bZnaedSzCS1NraqpycHKcPAcAy\nnVfuixcvJtQd4HiwL1q0SB9++KEk6cMPP1RFRYXThwBgoezsbO3bt49Qd0BSSzFLly7VgQMHdOHC\nBeXm5mr9+vVavHixlixZohMnTigcDuvTTz/VyJEjux+UpRgAGDTaHQHAMrQ7AvAEnS/eI9gBOIbO\nl/RAsANwBJ0v6YNgB5A0Ol/SC8EOIGl0vqQXdsUASBqdL6nBrhgAKdN552hJSUm35wl1bxDsABzR\nM9wJde+wFAPAUW1tbYrFYqquribUHcadpwBgGdbYASCgCHYAvVAL4G8EO4BuqAXwP4IdQAK1AHYg\n2AFIohbAJgQ7AEnUAtiE7Y4AJFEL4AdsdwQwKNQC2INgB5BALYAdWIoB0Au1AOmJSgEAsAxr7AAQ\nUAQ7AFiGYAcsRd9LcBHsgIXoewk2gh2wDH0vINgBi9D3AolgB6xC3wskF/exh8Nh3XbbbcrIyFBW\nVpYaGhr+f1D2sQOuoO/Fbp7foDRu3Dj98MMPGjVq1JCHAzB41wt3Qt0OaXGDEuENpB59L8h064ND\noZDmzJmjjIwMrVy5UitWrOj296qqqsTv0WhU0WjUrVGAwOkMd/pe/C0ejysejw/6fa4txbS2tmrM\nmDE6f/68SktLtW3bNs2ePfu/g7IUAwCD5vlSzJgxYyRJo0eP1iOPPNLty1MAgHtcCfbLly/rzz//\nlCT9/fffqq+vV3FxsRuHAgD04Eqwnzt3TrNnz9a0adM0c+ZMLVy4UGVlZW4cCrAenS8YLPrYgTTW\ndesiu1vg+Ro7gOTQ+YKhItiBNETnC5JBsANpiM4XJIM1diAN0fmC62GNHfCxnrUAnQh1DATBDqQp\nOl8wVCzFAGmura2NzhdISoPa3n4PSrADwKCxxg4AAUWwA4BlCHYgBeh7QSoR7IDLOvek19bWctco\nUoJgB1xE3wu8QLADLqHvBV4h2AGX0PcCr7CPHXAJfS9wGvvYAY/R9wKvEOyAi+h7gRdYigFSgL4X\nOIGuGACwDGvsABBQBDswQNQCwC8IdmAAqAWAnxDswA1QCwC/IdiBflALAD8i2IF+UAsAP2K7I9AP\nagGQTtjuCDiAWgD4kWvBXldXp8LCQk2cOFFvvPGGW4cBXEctAPzGlaWYjo4O3X333fryyy81duxY\nzZgxQx9//LEmT57830FZioEPUQsAr3m6FNPQ0KCCggKFw2FlZWXpySefVG1trRuHAlImOztb+/bt\nI9SR9jLd+NDTp08rPz8/8TgvL0/fffddt9dUVVUlfo9Go4pGo26MAgC+FY/HFY/HB/0+V4I9FArd\n8DVdgx0A0FvPi95169YN6H2uLMWMHTtWJ0+eTDw+efKk8vLy3DgUMGh0vsB2rgR7SUmJmpqa1NLS\novb2du3Zs0eLFi1y41DAoND5giBwJdgzMzO1fft2zZ07V0VFRXriiScSO2IAr9D5gqDgzlMEAneQ\nwgbceQp0QecLgoQrdgQCV+ywAVfsQBd0viBICHYEBp0vCAqWYhA4dL7ArwaanQQ7APgEa+wAEFAE\nOwBYhmCHL9H3AvSNYIfv0PcC9I9gh6/Q9wLcGMEO3+jr7lHCHeiOYIdv0PcCDAz72OEb9L0g6NjH\nDuvQ9wIMDMEOX6HvBbgxlmLgS/S9IIjoigEAy7DGDgABRbADgGUIdniGvhfAHQQ7PEHfC+Aegh0p\nR98L4C6CHSlF3wvgPoIdKUXfC+A+9rEjpeh7AYaOfexIS/S9AO4j2JFy9L0A7nI82KuqqpSXl6dI\nJKJIJKK6ujqnDwELdIb74sWLCXXAYY6vsa9bt04jRozQ6tWr+z4oa+wAMGierrET2gDgnUw3PnTb\ntm3avXu3SkpKtGXLFo0cObLXa6qqqhK/R6NRRaNRN0aBy6jPBdwTj8cVj8cH/b4hLcWUlpbq7Nmz\nvZ7fsGGDZs2apdGjR0uSXn31VbW2tmrXrl3dD8pSjBW6bl3kC1DAfWnRx97S0qLy8nIdO3ZsSMMh\nfV1vPzrhDrjLszX21tbWxO81NTUqLi52+hDwGLUAQHpz/Ir96aef1tGjRxUKhTRu3Djt3LlTubm5\n3Q/KFbuvVVRUqLa2ts+/L168WPv27UvhREAwpMVSTJ8HJdh9jVoAwBtUCsA11AIA6Y1gx5BQCwCk\nL5ZikBT2sQOpwxo7AFiGNXYACCiCHQAsQ7BD0n9r5RUVFdxcBFiAYEdiX3ptbS13jgIWINgDrufN\nRtQCAP5HsAcYnS+AnQj2AIvFYtetBZD+C/dYLJbiiQA4gX3sAUbnC+Av7GPHDdH5AtiJYA84Ol8A\n+7AUA0l0vgB+QFcMAFiGNXYACCiCHQAsQ7BbhL4XABLBbg36XgB0ItgtQN8LgK4Idp+j7wVATwS7\nz9H3AqAn9rH7HH0vQHCwjz0g6HsB0BPBbgH6XgB0xVKMReh7AexGVwwAWIY1dgAIqCEH+969ezVl\nyhRlZGToyJEj3f62adMmTZw4UYWFhaqvr096SADAwA052IuLi1VTU6MHHnig2/PHjx/Xnj17dPz4\ncdXV1WnVqlW6du1a0oMGBX0vAJI15GAvLCzUpEmTej1fW1urpUuXKisrS+FwWAUFBWpoaEhqyKCg\n7wWAEzKd/sAzZ85o1qxZicd5eXk6ffp0r9dVVVUlfo9Go4pGo06P4it99b2wbREIrng8rng8Puj3\n9RvspaWlOnv2bK/nN27cqPLy8gEfJBQK9Xqua7AH3Y36Xgh3IJh6XvSuW7duQO/rN9i/+OKLQQ8y\nduxYnTx5MvH41KlTGjt27KA/J0gG0veyb9++FE8FwK8c2e7YdV/lokWL9Mknn6i9vV3Nzc1qamrS\nvffe68RhrFVdXd2rEqBTSUmJqqurUzwRAD8bcrDX1NQoPz9fhw4d0oIFCzR//nxJUlFRkZYsWaKi\noiLNnz9fO3bsuO5SDP6PvhcATuLO0zTSda2dUAfQE5UCPkXfC4C+EOwAYBm6YgAgoAh2l1ANAMAr\nBLsLqAYA4CWC3WF9VQMQ7gBShWB30I2qAQh3AKlAsDtoINUAAOA2tjs6qK8rdom7SAEkj+2OHqAa\nAEA6INgd1jPcCXUAqcZSjEuoBgDgNCoFAMAyrLEDQEAR7ABgGYK9H/S9APAjgr0P9L0A8CuC/Tro\newHgZwR7D/S9APA7gr0H+l4A+B372Hug7wVAumIf+xDR9wLA7wj266DvBYCfsRTTD/peAKQTumIA\nwDKssQNAQBHsAGAZ64Pdi76XeDyesmN5gfPzN5vPz+ZzG4whB/vevXs1ZcoUZWRk6MiRI4nnW1pa\ndPPNNysSiSgSiWjVqlWODDoUXvW92P4PF+fnbzafn83nNhiZQ31jcXGxampqtHLlyl5/KygoUGNj\nY1KDJauvvhe2LQKw3ZCv2AsLCzVp0iQnZ3EMfS8AAs0kKRqNmh9++CHxuLm52dx6661m2rRp5sEH\nHzQHDx7s9R5J/PDDDz/8DOFnIPpdiiktLdXZs2d7Pb9x40aVl5df9z133nmnTp48qezsbB05ckQV\nFRX68ccfNWLEiMRr2MMOAO7pN9i/+OKLQX/gsGHDNGzYMEnS9OnTNWHCBDU1NWn69OlDmxAAMCiO\nbHfsegV+4cIFdXR0SJJ+/fVXNTU1afz48U4cBgAwAEMO9pqaGuXn5+vQoUNasGCB5s+fL0k6cOCA\n7rnnHkUiET3++OPauXOnRo4c6djAAIAbSPbL02S8++67prCw0EyZMsW88sorXo7imrffftuEQiHz\n+++/ez2Ko15++WVTWFhopk6dah555BFz6dIlr0dyxP79+83dd99tCgoKzObNm70ex1EnTpww0WjU\nFBUVmSlTppitW7d6PZLjrl69aqZNm2YWLlzo9SiOa2trM48++qgpLCw0kydPNt9++22fr/Us2L/6\n6iszZ84c097ebowx5rfffvNqFNecOHHCzJ0714TDYeuCvb6+3nR0dBhjjFmzZo1Zs2aNxxMl7+rV\nq2bChAmmubnZtLe3m3vuucccP37c67Ec09raahobG40xxvz5559m0qRJVp2fMcZs2bLFLFu2zJSX\nl3s9iuOefvpps2vXLmOMMVeuXOn3YsqzSoH33ntPa9euVVZWliRp9OjRXo3imtWrV+vNN9/0egxX\nlJaW6qab/vvHZ+bMmTp16pTHEyWvoaFBBQUFCofDysrK0pNPPqna2lqvx3LMHXfcoWnTpkmShg8f\nrsmTJ+vMmTMeT+WcU6dO6fPPP9ezzz5r3c67P/74QwcPHlRlZaUkKTMzU7fffnufr/cs2JuamvT1\n119r1qxZikajff5/Rv2qtrZWeXl5mjp1qtejuO6DDz7Qww8/7PUYSTt9+rTy8/MTj/Py8nT69GkP\nJ3JPS0uLGhsbNXPmTK9HccyLL76ot956K3HBYZPm5maNHj1asVhM06dP14oVK3T58uU+Xz/kSoGB\n6Gsf/IYNG3T16lW1tbXp0KFD+v7777VkyRL9+uuvbo7juP7Ob9OmTaqvr08858criIHcx7BhwwYN\nGzZMy5YtS/V4jguFQl6PkBJ//fWXHnvsMW3dulXDhw/3ehxHfPbZZ8rJyVEkErGyL+bq1as6cuSI\ntm/frhkzZuiFF17Q5s2btX79+uu/ITWrQ73NmzfPxOPxxOMJEyaYCxcueDWOo44dO2ZycnJMOBw2\n4XDYZGZmmrvuusucO3fO69EcVV1dbe6//37zzz//eD2KI7799lszd+7cxOONGzda9wVqe3u7KSsr\nM++8847Xozhq7dq1Ji8vz4TDYXPHHXeYW265xSxfvtzrsRzT2tpqwuFw4vHBgwfNggUL+ny9Z8H+\n/vvvm9dee80YY8zPP/9s8vPzvRrFdTZ+ebp//35TVFRkzp8/7/Uojrly5YoZP368aW5uNv/++691\nX55eu3bNLF++3Lzwwgtej+KqeDxu5a6Y2bNnm59//tkYY8zrr7/e705CV5di+lNZWanKykoVFxdr\n2LBh2r17t1ejuM7G/8R//vnn1d7ertLSUknSfffdpx07dng8VXIyMzO1fft2zZ07Vx0dHXrmmWc0\nefJkr8dyzDfffKOPPvpIU6dOVSQSkSRt2rRJ8+bN83gy59n479y2bdv01FNPqb29XRMmTFB1dXWf\nr/Xk/3kKAHCPfV8fA0DAEewAYBmCHQAsQ7ADgGUIdgCwDMEOAJb5H7IJs713AyMKAAAAAElFTkSu\nQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106672750>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 2"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.arange(1, 7)\n",
      "print('With numpy.outer')\n",
      "print(np.outer(x, x))\n",
      "\n",
      "print('The longer way')\n",
      "res = np.array([x * i for i in x])\n",
      "print(res)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "With numpy.outer\n",
        "[[ 1  2  3  4  5  6]\n",
        " [ 2  4  6  8 10 12]\n",
        " [ 3  6  9 12 15 18]\n",
        " [ 4  8 12 16 20 24]\n",
        " [ 5 10 15 20 25 30]\n",
        " [ 6 12 18 24 30 36]]\n",
        "The longer way\n",
        "[[ 1  2  3  4  5  6]\n",
        " [ 2  4  6  8 10 12]\n",
        " [ 3  6  9 12 15 18]\n",
        " [ 4  8 12 16 20 24]\n",
        " [ 5 10 15 20 25 30]\n",
        " [ 6 12 18 24 30 36]]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 3"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.outer(np.arange(1, 6), np.ones(2)).T"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "array([[ 1.,  2.,  3.,  4.,  5.],\n",
        "       [ 1.,  2.,  3.,  4.,  5.]])"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 4"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bucky = np.genfromtxt('../Data/bucky.csv', delimiter=',')\n",
      "print(bucky.shape)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(60, 60)\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 5"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.random.randn(5, 5)\n",
      "\n",
      "# x[0] is first row of x\n",
      "x[0] = [1, 2, 3, 4]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "ValueError",
       "evalue": "cannot copy sequence with size 4 to array axis with dimension 5",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
        "\u001b[0;31mValueError\u001b[0m: cannot copy sequence with size 4 to array axis with dimension 5"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y = np.random.randn(5, 6)\n",
      "np.concatenate((x, y))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "ValueError",
       "evalue": "all the input array dimensions except for the concatenation axis must match exactly",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-12-4d383539fa51>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[0;31mValueError\u001b[0m: all the input array dimensions except for the concatenation axis must match exactly"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 6"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "h = .001\n",
      "x = np.arange(0, np.pi, h)\n",
      "approx = np.diff(np.sin(x ** 2)) / h\n",
      "actual = 2 * np.cos(x ** 2) * x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.max(approx - actual)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "ValueError",
       "evalue": "operands could not be broadcast together with shapes (3141) (3142) ",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-14-6d4ff1880ba4>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mapprox\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mactual\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (3141) (3142) "
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "err = approx - actual[:-1]\n",
      "print(err.max())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.00987294960704\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(x[:-1], approx, x[:-1], actual[:-1], x[:-1], err)\n",
      "plt.legend(('Approx', 'Actual', 'Error'), loc=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "<matplotlib.legend.Legend at 0x106c16590>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD9CAYAAABDaefJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdUFOffBfC7sPSioDSpKhJBFBcRNKCuBTuICpbYS/Iz\ndk1sMcbeNcYSjUYjEjX2CGJvGKyoIHYRFUEUG1KlLLvz/uErCYoIW3h2h+/nHI6wMztznYXL8OwU\nAcdxHAghhPCCFusAhBBClIdKnRBCeIRKnRBCeIRKnRBCeIRKnRBCeIRKnRBCeEThUs/IyEBwcDBc\nXV3h5uaGixcvKiMXIYQQOQgVXcC4cePQuXNn7NmzB0VFRcjNzVVGLkIIIXIQKHLyUWZmJkQiER4+\nfKjMTIQQQuSk0J76o0ePYGFhgSFDhiA+Ph5NmjTBypUrYWhoCAAQCARKCUkIIVWNvPvbCo2pFxUV\nITY2FiNHjkRsbCyMjIywaNGij4Jp6sfMmTOZZ6D87HNUxfyanJ0P+RWhUKnb2dnBzs4OTZs2BQAE\nBwcjNjZWoUCEEELkp1CpW1tbw97eHgkJCQCAEydOoEGDBkoJRgghpOIUPvpl9erV6NevHwoLC1G3\nbl1s3rxZGbnUglgsZh1BIZSfLU3Or8nZAc3PrwiFjn757MIFAoXHhwghpKpRpDsV3lOXh7m5Od68\necNi1bxlZmaG9PR01jEIIYwx2VOnPXjlo21KCH8o8vNM134hhBAeoVInhBAeoVInhBAeoVInhBAe\noVIvg1gshrm5OQoLC1lHIYSQcqFS/4SkpCTExMTA0tISERERKlmHVCpVyXIJIVUXlfonhIWFoV27\ndhgwYAC2bNlS/PjgwYMxYsQItG/fHqamphCLxUhOTi6erqWlhdWrV6Nu3bqwsLDA5MmTiw9NCg0N\nha+vLyZOnIiaNWti9uzZyMrKwsCBA2FpaQknJyfMnz8fHMchPT0d9vb2iIyMBADk5OTA2dkZW7du\nrdwNQQjRLJwKfWrxKl6tUtStW5fbunUrl5CQwOno6HAvXrzgOI7jBg0axJmYmHDR0dFcQUEBN27c\nOM7Pz6/4eQKBgGvTpg335s0bLjk5mXNxceE2btzIcRzHbd68mRMKhdyaNWs4qVTK5eXlcQMGDOCC\ngoK4nJwcLikpiXNxceE2bdrEcRzHHTt2jLO2tuZevHjBDR8+nAsJCflkXk3YpkSzpLzI4n4MO8D1\nX7GBm/PXQe51Vi7rSFWGIj/PalvqgOIf8oqOjub09fW5rKwsjuM4zsPDg1uxYgXHce9KvW/fvsXz\n5uTkcNra2tyTJ084jntX6kePHi2evnbtWq5t27Ycx70rdQcHh+JpRUVFnK6uLnfnzp3ix9avX8+J\nxeLir8eMGcO5u7tzdnZ2XHp6ehnbi0qdKEdmTgHX4sfZnGCqGWc+vi1X7/uhXLWxYk4wxZzruexn\nrkgqZR2R9xT5eVbb4Rdl1Lq8tmzZgvbt28PExAQAEBISUmIIxs7OrvhzIyMjmJub4+nTp8WP2dvb\nF3/u4ODwyWmvXr2CRCKBo6NjiflTU1OLv/76669x69YtDB48GGZmZvL/pwgphxsPXqLWtDZIeHsJ\n//S/gtcrTiBh6SZkrDyN8G7ncTR5L2pPCUHW23zWUcknqG2ps5KXl4ddu3bh1KlTsLGxgY2NDZYv\nX47r16/j+vXrAICUlJTi+XNycpCeno5atWoVP/bfMfbk5GTY2toWf/3fu0HVrFkTOjo6SEpKKjH/\n+18aUqkU33zzDQYOHIhff/0VDx48UPr/l5D3ElLeoOnq9mhcwxepSw/Ar0GdEtMDmn+BlHknISkU\noNHMASiiN/rVEpX6B/bv3w+hUIg7d+4gPj4e8fHxuHPnDvz8/BAWFgaBQIBDhw7h3LlzKCwsxIwZ\nM9C8efMSxb1s2TJkZGQgJSUFq1atQu/evUtdl7a2Nnr16oXp06cjJycHjx8/xooVK9C/f38AwIIF\nC6CtrY3Nmzdj0qRJGDhwIGQyWaVsB1K1vM0vQtOlPdHItBWif1oEba3Sq6G6iR5uz92G9IIX6LJk\nXiWnJOWixGGgj3xq8SperUI6duzIff/99x89vmvXLs7a2prr378/N2LECM7f358zNjbmWrVqxSUl\nJRXPJxAIuNWrV3N16tThatSowX3//fecTCbjOI7jQkNDuRYtWpRY7ps3b7j+/ftzFhYWnL29PTd3\n7lxOJpNxV65c4czMzLgHDx5wHMdxUqmU8/X15RYsWFBqbnXepkT9eU3+gbOY2I6TFBWVa/4LN1M5\nwWRLbsup8ypOVjUp8vNMV2msoCFDhsDOzg5z584tdbqWlhYSExNRp06dUqeriiZvU8LWyn3n8d2l\nYNweew0utpblft6INTsR9mg+MhbHQlfI5CrevEVXaaxEVJyETzKyCzE5+htMavhLhQodAH79theE\nkhr4+rcNKkpH5EGlXkECgaDEm52lTSdEU3Rf8jNqaDtiQb+QCj9XW1uAX7uuxNaU2XiVlaOCdEQe\nNPzCE7RNSUXF338J0SZXnB0Ugy9d5R8urDW2N5o7NsXe775XYrqqjYZfCCEV1m/9fIh0+ipU6ACw\nuOt0hL9Yjqy3eUpKRhRBpU5IFXTyahJuC//E9hE/KrysAe0bofpbb4zfHKaEZERRVOqEVEHfbluM\nlkb/wxe2VkpZ3ljvsdj5cC0NAaoBKnVCqpgrd58jUW8Hfh82XmnLnNanDQpl+dhy6rzSlknko3Cp\nS6VSiEQiBAQEKCMPIUTFRm1ZjUbafVCvVsUOYSyLjo4A/mbfYt6xX5W2TCIfhUt95cqVcHNzo0P5\nlGzWrFkYMGAA6xiEZx6nZeMy9xvWfPWd0pe9fOAgPNA+iKfpGUpfNik/hU4De/LkCQ4dOoTp06fj\n559/LnWeWbNmFX8uFoshFosVWWWlEovFuH79OtLS0qCrq1vmvKGhodi0aROio6OVsm76JUlUYdTv\nm+Egaw0/N2elL9vVyQyWOe0wa9cebBgxXOnL57OoqChERUUpZVkKlfqECROwdOlSZGVlfXKe/5a6\nJnl/OzsHBwdEREQgODi4UtdPbzgRZSsq4nD09W9Y12WdytbRx60/tiWuxAZQqVfEhzu8s2fPlntZ\ncg+/REZGwtLSEiKRiJcF9Knb2aWkpKBHjx6wtLREzZo1MWbMGNy9excjRozAhQsXYGJiAnNzcwDv\nXqhNmzYVPzc0NBQtWrQo/nrcuHFwcHBAtWrV4OXlhbNnz1bef5BUOUt3/wMdoQDD2rVU2Tpm9O6M\ndOENXHuU/PmZiUrIXernz59HREQEateujb59++LUqVMYOHCgMrMxFRYWht69e6NXr144evQoXr58\nCalUiq5du6J27dp4/PgxUlNT0bdvX9SvXx/r169H8+bNkZ2djfT0dACfv6SAt7c34uPj8ebNG3z1\n1VcICQlBYWFhZf0XSRWz+sI69HAYodKhvZpmeqhT0BOz9v6lsnWQssk9/LJgwQIsWLAAAHDmzBks\nW7YMYWHKO/lAMFvxbzxupnx/QZw9exapqakIDAyEiYkJ3NzcsG3bNvj4+ODZs2dYunQptP7/etNf\nfvnlu3XJ8ddKv379ij+fOHEi5s2bh3v37qFhw4Zy5SbkUy7feY40o6NY0v83la+rvygEK+JnAJii\n8nWRjyntepnK/u0vbyErw6duZ2drawtHR8fiQlfUsmXL8Mcff+Dp06cQCATIysrCq1evlLJsQv5r\n0vY/4CboiVrm1VW+rgndxZh9OwG3UlLRwN72808gSqWUUm/VqhVatWqljEUx9/52djKZDDY2NgCA\ngoICZGZmwsrKCsnJyZBKpdDW1i7xvNJ+qRkZGSE3N7f467S0tOLPo6OjsXTpUpw6dQoNGjQAAJib\nm/Py/QnCVlERh7O5f+DPoK2Vsr5qJjqwy+uCpRHhCB01slLWSf5FZ5R+oKzb2f3999+wsbHB1KlT\n8fbtW+Tn5+P8+Xdn0FlZWeHJkyeQSCTFy2rcuDH27duHvLw8JCYmYtOmTcXln52dDaFQiJo1a6Kw\nsBBz5swp8ygiQuT1a8RFCLW00aeFd6WtM+iL7jj86O9KWx/5F5X6B8LCwjB06FDY2dnB0tISlpaW\nsLKywujRo7Fz505ERkYiMTERDg4OsLe3x65duwAAbdu2RYMGDWBtbQ1Ly3dn6k2YMAG6urqwsrLC\nkCFDiu89CgAdO3ZEx44d4eLiAicnJxgYGMDBwaF4+ufeZCWkvH49G4Z2FgMr9ftpUo8OeKEbg6dv\n3lTaOsk7dD11nqBtSkrz6k0BLBfXwtVv4iCq4/D5JyhRzTEBGNa0PxYPLP3G6+TT6HrqhJBSzf4r\nEjWKPCq90AGghU0nRNw5Uunrreqo1AnhsR13wtC7PpvzR4aLO+C+7Aj9BVnJqNQJ4alrCa/w2uQM\nZvXqyWT9nZvXBSTGOHj1OpP1V1VU6oTw1Ny9e1BH2hk1TU2YrF8gAOqhIzZG0RBMZaJSJ4Snjj/d\nhQGevZhmCHTriLNpVOqViUqdEB6KuZ2GHOM4fNetI9Mco7qK8VrvCl5nZzPNUZVQqRPCQwv374Mz\n1wXG+vpMczhYG8EkywsbjyvnPgPk86jUCeGhE892YZAX26GX9xqZtkbEjdOsY1QZVOqE8MyFm8+Q\naxKP8QHtWUcBAAQ2EuNGVhTrGFUGlfoHnJycYGhoCBMTk+KPsWPHso5FSLkt3L8X9WQBMNJjO/Ty\n3tD2PsjWv4PnmZmso1QJVOofEAgEiIyMRHZ2dvHHqlWrPppPKpV+9JhMJqvQuio6PyHlcer5Lgxu\nGsI6RrGaZnowzfLBHydoXL0yUKmXU2hoKHx9fTFx4kTUrFkTs2bNwpAhQ/Dtt9+ic+fOMDY2RlRU\nFO7cuQOxWAwzMzO4u7vjwIEDxcsYPHjwR/MTokzR157irckNjOuqHkMv73lUa40DNK5eKajUS/Gp\n05pjYmJQt25dvHjxAtOnTwfHcfjrr78wY8YM5OTkoGnTpggICEDHjh3x8uVLrF69Gv369UNCQkLx\nMv47v6+vb2X9l0gVsSRyL1y4ABjq6bGOUkI3DzFu5ESxjlElqG+pCwSKf8iB4zgEBQXBzMys+GPj\nxo0AgFq1amHUqFHQ0tKCvr4+BAIBgoKC0Lx5cwDAtWvXkJubi6lTp0IoFKJ169bo2rUr/vrr3/s1\n/nd+PTX7wSOa758Xf+OrxmwuC1CWwf7eyNFLwLMMuhSvqqlvqXOc4h9yEAgECA8Px5s3b4o/hg8f\nDgCwt7f/aH47O7viz58+ffrRPI6Ojnj69GnxsktbBiHKcPtROrKMr2BcgD/rKB+pUV33/8fVz7KO\nwnvqW+pqqLSbDPz3sVq1aiElJaXE8M3jx49ha0v3aSSqtyz8EGwLW6OaoSHrKKVqYOqLo7fPs47B\ne1TqpSjvpUI/nK9Zs2YwNDTEkiVLIJFIEBUVhcjISPTp06dCyyVEHocfhqNrvW6sY3xS+/q+uJlJ\npa5qVOqlCAgIKHGceo8ePUq9vdyHj+no6ODAgQM4fPgwLCwsMHr0aPz5559wcXEpdX5ClOXlm3yk\nGR3HpG5dWUf5pP7iZnijfxX5kkLWUXiNbmfHE7RNq7YfNh/CupsL8Ga5eo9Z643zQGjPDejb0od1\nFLXG7HZ2KSkpaN26NRo0aAB3d/dST9IhhKjenhvhaFMriHWMz3LS9sX+qzQEo0oKlbqOjg5WrFiB\nW7du4eLFi/j1119x584dZWUjhJRDQaEMicIITOikvuPp7/na++LSs3OsY/CaQqVubW2Nxo0bAwCM\njY3h6upafPgeIaRybDp8GXoyM/i51WMd5bNCmn+JJ1rnaKhQhYTKWlBSUhLi4uLg41NyrGzWrFnF\nn4vF4uJT6OkNQ+UyMzNjHYEwsvlCOJqaqv9eOgD4ezlBtleAqw+T4FW3Nus4aiMqKkpplw1Ryhul\nOTk5EIvF+PHHHxEU9O+4Hr15R4hqcRygN6EBNnffhH6tmrGOUy5WY4LRr0kQfh7cn3UUtcXsjVIA\nkEgk6NmzJ/r371+i0Akhqhd5PhEyvXT0aeHNOkq5iWr64nQivVmqKgqVOsdxGDZsGNzc3DB+/Hhl\nZSKElNOvJ8Lhqh0AbS3NOeWko3szJObFsI7BWwp9J5w7dw5bt27F6dOnIRKJIBKJcOQI3TmckMpy\n7nU4+jXRjPH093q1bIwcg9vIyc9nHYWXFHqj1M/Pj270QAgj8fdfIdc4HqO7tGUdpUJqWRhAL7s+\nwi9d05j3ATSJ5vzNRggpYVlEJOwl7WCsrx63rasIey1vRMbREIwqUKkToqGOPg5Ht/qaNfTynpeN\nN66mUamrApU6IRoo7XUeXhqfwneBXVhHkUuXxt5IllKpqwKVOiEaaPnfJ2BeIIKjRQ3WUeQS5OuK\nAp00PH1Dd0JSNip1QjTQvlv70c5OM4deAMDYSBvG2SLsPneFdRTeoVInRMPk5UvxSCcSE7tobqkD\nQF19bxy/RUMwyqa0a7+QzyuSynA14SniH6XgwfM0PH6dhqfZz5BRkI58aR4KZG9RKHsLGWQQCnQh\nhC6EWrqopmsOC0ML2Jhaop6VHVq4uaBZfUfo6miz/i8RBn6LvAgDmRV8XOqwjqKQ5g7eiEzeyjoG\n71Cpq4CkSIbjV+/jYFwsYp9cR3LufbxGAgoMH0CryBSGhY4wEVjDTMcalobWqGfuAhN9IxjrGcBE\n3xDaWlookEiQLylEnqQAL7LT8SL3BR5lPEDEw2TMvJwAqf4L6L2tC0ehN5rbNUcP7+bo4u2mUWcW\nEvmExYTDp7pm76UDQJC3N35PHQOO4+gCf0rE5M5HfMJxwLmbydh+9izOJl1EUkEsso3iISywgBUn\nwhfVPNCo1hfwdq4HcaN6sDE3Ucp6X2a8xZErdxF57RIuP7uAFJyHTJiLelwXdHfvgkndO8DcRD1v\nQEzkx3GA7sQv8FfINgR/6cU6jkIkEg66061xY8wVuNvbs46jVhTpTir1CpJKOYRfuI3dl6JxIfUs\nnmhFgxPmwUbSAp4WzdG6vid6NBfB0apyL4XLccCxK4lYd+Ig/kmLRIbhFbjIumN0i4H4tnNL2oPn\nib1n7qLPoXYoXJTCi73bGmO6YozfUMzq3YN1FLVCpa5i1xKfY+2RYzj28BiShcehLTOEk6Al/Bxa\noK+vH9o1doGWlnr9gF1NeIafdm/HyVeh4LQk6OM4ASuHDkR1YwPW0YgC2s5cjHTpY8TNW8s6ilI0\n/2Em9A2kOD1jHusoaoVKXckysguw/shZ7Ik7hpt5R1Gg/xi1CtugjUN7fNOuPfwaaM7F/WUyDj/v\nO4Ml0T/jtcFF+JuOxZ+jxsOimjHraEQORmO/xNx2MzExsAPrKErx3e/h2HZ3PdKWH2IdRa1QqSvB\n7aRXWB5xEIcfRuCZ4QmY5DWAl1l79PVuj4FtvaGno/nvKR+4cA8jd87GU73TCLH6ARu//QbGBnqs\nY5FyirmdhmZb6yNn5nMY6vHjdYuKS0HbnV4oWpjGi+EkZaFSlwPHAYdjErDmeATOvoxAtlE8bAva\noatzIL4L7IJ6tjVZR1SZbSevYVzED8jWScRCv7WYGNSOdSRSDr2X/I7Lr07i4ZIdrKMojVTKQecH\nK9wYHYcG9ras46gNRbpT83c/K6CgUIrfj1zElosRiM8Ph1SYBVetQEz5chrGBLSGqaHmXe1OHv3a\nNsZXbQ7hh9BITD77Ndaca4aIUT/D3cmGdTRShpNPwjGkST/WMZRKW1uAam89sf9iLJW6kvB+Tz2v\noAirI/5BaMxu3BX8Db0iS3hX64ZvWnZDn1aeVf6okOfpbxGwdC6uyv7AJPc1WDQghHUkUorHz3Lg\ntLoWnk5Kho1ZddZxlKrp1B9gbqqHoz/MZB1FbdCe+gfe5hdhxf7T+PPKHiRo/w1DiQNaWYRgbZez\nEDdyZh1PrViZGyJm4UKsjwzC6FMDsO/mAZyevBq2Naqxjkb+Y/n+Y7Ao9OFdoQOAj4Mn9j8KYx2D\nN3hT6jlvJVix/xT+jN2NRO1wGBXWQWurYGwMuKRRR6uw8r+uPghsHgfxgu9Qe6EndgXvQ1AzD9ax\nyP8Lv7cfHevy88buXTybYP1jusexsmj08EtWbiGW7TuB7df24KFOOIwLXNDWJgRTAnuiWX1Hla2X\n775Z/Rc2po7FWJdf8MtQfo3haqLs3CJUm2uNayPi0MiJf2deSiQcdH8yR+K4u6hrbcU6jlqoUsMv\nGdkFWPb3cWy/thtJegdgku+GdrWCsaPbbHi58O8bnoUNY/qi9akGGHiwB87/dAXnfloGHSFdPIyV\nNRFnYVTkyMtCBwAdHQFMczyx/1IcvuvWkXUcjacRe+rpWflYvPcodt7Yg8d6kaiW3xAd7EIwNagH\nRHXpHXNVuZ/yBk2XhsBY1wg3Zv0FM2O6lgwL7t+Ph33NGjg8dQbrKCrTePIk2NWojsgp01lHUQu8\n3FN/lZGHhXsOY/etPUjRPwSzfBE62AfjQPclaFibDr2rDPXszfB4wSF4/Pg1HH9qjcsTD+ALO0vW\nsaqUoiIOd7hwLPEPZx1FpbztmuDok92sY/CCWpX6izdvsXDPIey+vRup+kdhnu+Fzk7BmNZ9Bdwc\naKyNhWrGuniwPBQtZsxCw1++xD/DT6BZfSfWsaqMbSduQFsb6OTZkHUUlero4YnQJ9NYx+AFhQ/S\nPnLkCOrXr4969eph8eLFFX7+01c5GLt+F2wnhMBqqQ3+vL0B7eq0xd1R9/F6xQn8OW4EFTpj2toC\nnF8wG53MJqDFplaIunGfdaQq4/focHgaBPH+FPpOPs6QCF8jNT2ddRSNp9CeulQqxejRo3HixAnY\n2tqiadOmCAwMhKura5nPS3mRjYV7I/H3vd1IMzwJi/wvEVA3GNN6rINzLf6enq/pwqePQp8lemj3\nZ2sc7HMcHTzLfp2JYjgOuJK9H6u6LmcdReUM9LVglC3C/otxGNW5Les4Gk2hUo+JiYGzszOcnJwA\nAH369EF4eHippf44LRPz9xxAROIePDc8Bcv8FujmEoIfem6Ek5W5IjFIJdoxeTiGrNBD5x1tESE9\nhi5N3VlH4q3TsSmQGD3GkLZ+rKNUCic9T5y8c5VKXUEKlXpqairs/3PHEjs7O1y6dKnEPO29/fAk\n+zHyhWlwMHfBVFFX9G0yBFbVjN7NcCMWuPHBgst611fTp6lLDgWmbbY3gPf1XtgyryWcWs1GAweb\ncj2PtmXFXI84iUkFX0Bn67aPJ8ozHCMQANragJaW8v/98DGB4N3npf37iWnNLNxx7ukRQCL5eH6O\nA2Sydx9Safn/ff8hkXz6o6io5Ne9egGVfBXMqKgoREVFKWVZCpV6ecb5ZqTdg6WRJRzN6kBfKARu\nXHn38W4BZS2cv9PUJYcC074VCCB6WBt3705FrTpimBkZlOt5tC3Lz/L2cYgt6wOnTpWcIO9hwvKU\nYkX+/W+Rcty/RVzav6U8tl4qhVQqAVbtLznff7ehPL90tLUBHZ2yP4TCfz/v1q3SS10sFkMsFhd/\nPXv2bLmXpVCp29raIiUlpfjrlJQU2NnZlZinRfJLRVZB1FgzAB1n/4x+2RtwfcI/cLGlwx2VJfbe\nS/TfUg8ZPx4CDKvG3apyc4pQbWE1vJqahhom/7mXr0z27x47+SyFjn7x8vLC/fv3kZSUhMLCQuzc\nuROBgYHKykY0wOGfJsJTtzdEP3fAi8xs1nF4Y9H+/agt7QjTKlLoAGBqLIRBtjsiYuJLTng/DEPK\nRaFSFwqFWLNmDTp06AA3Nzf07t37s0e+EH4RCIDoubNgI2uKRnN7oUAiYR2JF44/2YM+jXqyjlHp\n7LQ9ceJWLOsYGk0jLhNA1F9uXhHspwTAzsQe8fPW8/64alW6+eA1Gm6qjdc/PIW5cdW6l2zIog24\nlXkBtxduZh2FKUW6s2rfIYIojZGBENem70JC7mUELl3EOo5GW/h3BByK/KtcoQNAG1dPJEtoT10R\nVOpEaRysTHDmm4M4/OI3TNrCn/toVrYjj/eil3sw6xhMBH3pjlz9+3hbmM86isaiUidK5eNWC1s6\nHsDyO2Ow5xztcVVUwuNMpJtEY3JQF9ZRmLCx0IdOdj0cvnqTdRSNRaVOlK5fu0YY6bAOffd3R0Lq\nC9ZxNMrCfQdgK2kFC1NT1lGYsYEIR6/HsY6hsajUiUqsGRkMT+FA+CwPRl5hIes4GiPy4R70dK2a\nQy/vudfwxOUU+itPXlTqRGWi58yGsKg6ms8ZxzqKRniYmoVXJqcxpXsA6yhMtXIR4VE+7anLi0qd\nqIyujhYuT9uK27lnMHzdBtZx1N7snfthWyhGLTMz1lGYCmrWGJl6NyGRFrGOopGo1IlKOdmYYk/w\nfvzxeDp2n7vKOo5ai3i0Hf0afcU6BnP1HE2g/bYW/rl1j3UUjUSlTlQu0NcFI+zXot/+EKS8esM6\njlq6evcFMk0uYlrPqj30Arw7S9miSISDsTQEIw8qdVIp1o4KQT1ZAJovHgyZjM4y/tDcvbvhLOuK\n6kZ0c28AqF/dExcf05ul8qBSJ5Xm3KyleFP4HMErlrGOonaOP9+O4T409PKebx0R7mfTnro8qNRJ\npaluootDQ3Zh//Pl2HgsmnUctXH0UhLyDRMwrqs/6yhqI9BbhNe6cXTtKDlQqZNK1aqxA6a5bsaI\nE33pxKT/tyhyBxoKg6Gno8M6itrwcrUACk0Q++gR6ygah0qdVLr5QzqhsWAgWq4YCKlMxjoOUxwH\nnMvajrGtaejlv7S0APN8T4TH0Lh6RVGpEyaiZs5GdmEWev/yM+soTP1xOA7Qy8Kg1r6so6gdZ2MR\nzj6gcfWKolInTBgb6uDA4O3Y93wJdvxzmXUcZpaf/APi6oOhrUU/ih9q5uiJuxm0p15R9J1EmGnj\n6YQRDr9i0IG+SHuTxTpOpXudUYC7On9hfshg1lHUUhdPEV5o0556RVGpE6Z+HRkC+6I2aLl4ZJU7\n0mHGtgjUkHigaT0n1lHUUiuRHaScFIlpz1hH0ShU6oQpgQA4++MvSCqIw9hNYazjVKoddzejv/sQ\n1jHUlq6mNSkRAAAStUlEQVSuAKa5Ivx9iYZgKoJKnTBnXcMQoV134NfE73H6egLrOJUi+loqMowv\nYlbvHqyjqLXa+p44c4+GYCqCSp2oha/aNkQ30zkI2NIHufkFrOOo3Iy9YXDXCkY1Q7osQFm87ES4\n+Zr21CuCSp2ojT1TRsBI4oS2C6eyjqJS+QVSnH37O6Z3Gs46itrr0EiEZxztqVcElTpRG9raAkRN\n3IgrufuwcM9B1nFUZtbWIzBEDfT282YdRe11buaMQuFrpGWms46iMeQu9UmTJsHV1RUeHh7o0aMH\nMjMzlZmLVFGuTuZY0mwbfrw8DNcfpbKOoxK/X1uL/vVHso6hEYwMtWCY5YHwS9dYR9EYcpd6+/bt\ncevWLcTHx8PFxQULFy5UZi5ShU0M9oOvzmi0XtMPkiIp6zhKdfjiQ2QYxWBRvz6so2gMex1PnLxN\n4+rlJXep+/v7Q+v/z4Lz8fHBkydPlBaKkGMzpqFIoo2uS+exjqJUP+z7DT76g2BqaMA6isYQWYsQ\n/5zG1ctLqIyF/PHHH+jbt2+p02bNmlX8uVgshlgsVsYqCc/p62nj+IitaB7qid+OiDGiYyvWkRT2\nPD0P8YLNiO59gXUUjdKugSfCzyxmHUOloqKiEBUVpZRlCbgyTuPz9/dHWlraR48vWLAAAQHvbrs1\nf/58xMbGYu/evR8vXCCocmcJEuX6MfQIFt3+GvcmxKKujQXrOAoJWbQe514dwNNlkayjaJSXryWw\nXFENmdNfwtTAiHWcSqFId5ZZ6p8TGhqK33//HSdPnoS+vr5SgxHynmjyFDyX3cSTJZHQ0hKwjiOX\n/AIpTKa5Ym3Hjfi6fUvWcTSO/pimCO29Cn38mrOOUikU6U65x9SPHDmCpUuXIjw8vNRCJ0RZ/pk5\nD5mS1+j9ywrWUeQ2ZXMEDAVmGO7fgnUUjVRLS4SjN+jN0vKQu9THjBmDnJwc+Pv7QyQSYeRIOkSL\nqIaJkQ4iBu7A3ueLNPIyvTIZh99vL8HYJpMhEGjmXxqsNazpiaup9GZpecj9Run9+/eVmYOQMrVt\n4oSRF9dh4IE+aNUwFjZm1VhHKrfle85CqvcKM3sHsY6iscRfiHD68gbWMTQCnVFKNMbqkT1RW9oB\nvgu/0Zj3ajgOmBc9G4PqToFQW5t1HI0V9GUjZOvdRUFRIesoao9KnWgMgQCI/nE5nhUkos8vv7CO\nUy6Ld55Bnv4jrB42iHUUjVbbzgDC7Do4eeMW6yhqj0qdaBRLcwMcHLQPu58txrrDUazjlEkm4zD/\n3E/4X/2foKejwzqOxrOUiXA4jsbVP4dKnWicNp6OmNnoT4yO6ovYByms43zSnK0nIdFLw8+D+7GO\nwgtuZp64lExHwHwOlTrRSDP7+0OsNx7itT2Rk5/POs5H8vKlWHR1EiY0ngsdbaWcuF3ltXAWITGX\n9tQ/h0qdaKyjP02GocQBzeaMVrs3Tgev3AwDoTHmfxXCOgpvBHo3RoZePKQyfl3kTdmo1InGEgoF\niJm+GfffxiD45+Ws4xR7mJqF3a9nYEOPXzT2DFh15PFFdQhyrXApkQ6nLguVOtFoDlYmODn0IPan\n/YLp2z6+/hAL3VbMQn3tTgjxbcI6Cq8IBECNQk8cuEJDMGWhUicaz6+RPX5vE4GFN0Zg6+lLTLOs\nDb+E29rbETluCdMcfOViKsL5h/RmaVmo1AkvDO3kiUn1NmPQ4e44cyORSYb0zAKMjxqK7xr8gjrW\nNZlk4LvmTp64m0l76mWhUie8sXhYV4TUnIN2Ye1wOeFxpa/ff8FMWAqdsXhA70pfd1XR1UuEVzqx\navfGuDqhUie8smPycHQwnQjfDW1xI+lppa13+ubDiOe2IWrCRrpolwr5eliBKzTArdRk1lHUFpU6\n4Z3IGWPhazAMTde0RdwD1d9mMTo+BQvvDMHattvgXEuzb+Sh7oRCoFqeCOExNK7+KVTqhJdOzZmG\nliZD4f2bH07FJ6hsPQ9TM9EutDN6WE/CNx3o5heVoY6BJ6Lv07j6p1CpE14SCIBjMyehl9VP8N/e\nCn+q4KiY15n58FwcjAbGLbF7wkSlL5+UzttBhFvptKf+KVTqhNe2fT8U33+xAYOOBGD42o1KW27a\n67dwnhEIc/0auDhzJY2jV6KOHiKkCWhP/VOo1AnvLR4agPBu0Qi7/zPqTx2C1NeZCi3vyt1ncJ7T\nDpaGNri3YBt0hXRtl8rU3tsRRchDcvpz1lHUEpU6qRICvvwCD6ZdAiT6cFzsjp+274dMVvHD4uZu\nPw6fTU3ha9URtxdsho6QbnxR2QwMBDDKFmH/RdpbLw2VOqky7C1NcHf5OsxrEobFl6fD7Hs/LNl7\nDFKZ7LPPjbx4D/bjv8KcuK+xqPkmHP3hJ2hr0Y8PK3UNPHHiFpV6aQScCo/iFwgEdJIAUUsFhVKM\nXv8Xtj5chiJhBkR6wejo2gp+9b+As40FsvMKcDUxGYfjr+Bk6n5kGlxDG+NR2DH2e9QwNWIdv8r7\n3+q/cDh5D5KXqsf1fpRNke6kUidVGscBf56IRej5g4h/E40s7Yco0n0FgUwP+hIb2AlFCHLrgik9\nOlGZq5Fjlx+h8x4/SBY94eWb1FTqhJAqRSrloDPNBldHxEBUx4F1HKVTpDtpUJAQonG0tQWomd8M\nO85eYB1F7Shc6suXL4eWlhbS09OVkYcQQsqlYfXmOJ1Ipf4hhUo9JSUFx48fh6Ojo7LyEEJIubR3\na457uVTqH1Ko1CdOnIglS+hmAISQyte3lReyDG4it0D9bjzOktynwoWHh8POzg6NGjUqc75Zs2YV\nfy4WiyEWi+VdJSGEFHOwMYReliv2XbiKAWJf1nEUEhUVhaioKKUsq8yjX/z9/ZGWlvbR4/Pnz8eC\nBQtw7NgxmJqaonbt2rhy5Qpq1KhRcuF09AshRIXqTxiDhg6O2D3he9ZRlEqR7ixzT/348eOlPn7z\n5k08evQIHh4eAIAnT56gSZMmiImJgaWlpVxBCCGkoprZNceZNH6egCQvpRynXrt2bVy9ehXm5uYl\nF0576oQQFTp0/hECw30hWZTKq5OQmB+nzqeNSQjRHO29nSDjZLiaSLe3e08ppf7w4cOP9tIJIUTV\nhEIBrAtaIjTqDOsoaoPOKCWEaDQfKzFOPaRSf49KnRCi0fo0E+NBURTrGGqDSp0QotG6+7lCopWN\n649pXB2gUieEaDhdXQEs3rbC5tM0BANQqRNCeMDbQowT96nUASp1QggPhDQV474kinUMtUClTgjR\neL3buKFQkIVrD1NYR2GOSp0QovH09ASwKWiN346dYB2FOSp1QggvtLbviKMPD7OOwRyVOiGEF771\n74jH2icgkRaxjsIUlTohhBe+bGgDnbcO2BF9iXUUpqjUCSG8IBAAbrqdsOV81R6CoVInhPBGL1En\nxKQfYR2DKSp1QghvjApsjmydB7iR9JR1FGao1AkhvGFqrAOHvK5YHP436yjMUKkTQngl2C0ER1J2\ns47BDJU6IYRXJvVoj9c68UhMS2MdhQkqdUIIr1jX1IdNTmfM31s1h2Co1AkhvNPPoxf2P9jOOgYT\nVOqEEN75qW9nZArv4/SNu6yjVDoqdUII75gY6aARNwg/7d/EOkqlo1InhPDSzIDhOJ8bhpy8gnI/\nJ3DWJiQ/z1ZhKtWjUieE8FL3lvVQ7W1jTNyytVzz7zx9CwfzfkTN6roqTqZaCpX66tWr4erqCnd3\nd0yZMkVZmQghRCkmNfsBYQ8WoUgq/ey8P+xfhfZm38JQT68SkqmOUN4nnj59GhEREbh+/Tp0dHTw\n8uVLZeYihBCFTe7dEvMuWmLC5u1YPXzAJ+c7FZuERwZ7cHKo5r+xKvee+rp16zBt2jTo6OgAACws\nLJQWihBClEFbW4BFrZdh3f2peJGZ9cn5BoXOQGvj0XCy1PweE3Acx8nzRJFIhG7duuHIkSPQ19fH\nsmXL4OXlVXLhAgFmzpxZ/LVYLIZYLFYoMCGEVFSd8cNhaAjcXLDxo2kzthzGopvf4umM67AwNWWQ\nDoiKikJUVFTx17Nnz4ac1Vx2qfv7+yOtlFNt58+fj+nTp6NNmzZYuXIlLl++jN69e+Phw4clFy4Q\nyB2MEEKUJfl5NpwX+aCH4/+wY/y44sfDz95DjwgxVou3Y2Tn1gwTlqRId8q9p96pUydMnToVrVq1\nAgA4Ozvj0qVLqFGjhlKCEUKIMp2Ke4QOf7ZHHR0/DPHqi0uJ9xCeMQ8j6y3BmuGDWMcrQZHulHtM\nPSgoCKdOnQIAJCQkoLCwsEShE0KIOmkjqo3EyVdgqeuIxecWIf71JWztHKl2ha4ouffUJRIJhg4d\nimvXrkFXVxfLly//aLyc9tQJIaTimAy/lGvhVOqEEFJhTIZfCCGEqB8qdUII4REqdUII4REqdUII\n4REqdUII4REqdUII4REqdUII4REqdUII4REqdUII4REqdUII4REqdUII4REqdUII4REqdUII4REq\ndUII4REqdUII4REqdUII4REqdUII4REqdUII4REqdUII4REqdUII4REqdUII4REqdUII4REq9TJE\nRUWxjqAQys+WJufX5OyA5udXhNylHhMTA29vb4hEIjRt2hSXL19WZi61oOnfGJSfLU3Or8nZAc3P\nrwi5S33y5MmYO3cu4uLiMGfOHEyePFmZuQghhMhB7lK3sbFBZmYmACAjIwO2trZKC0UIIUQ+Ao7j\nOHme+PjxY/j5+UEgEEAmk+HChQuwt7cvuXCBQCkhCSGkqpGzmiEsa6K/vz/S0tI+enz+/PlYtWoV\nVq1ahe7du2P37t0YOnQojh8/rpRQhBBC5CP3nrqpqSmysrIAvCvv6tWrFw/HEEIIYUPuMXVnZ2ec\nOXMGAHDq1Cm4uLgoLRQhhBD5lDn8UpYNGzZg1KhRKCgogIGBATZs2KDMXIQQQuQg9566l5cXLl26\nhGvXrmHmzJno27cv6tWrh8WLF5c6/9ixY1GvXj14eHggLi5O7sCqcOTIEdSvX/+T+aOiolCtWjWI\nRCKIRCLMmzePQcrSDR06FFZWVmjYsOEn51Hnbf+5/Oq87QEgJSUFrVu3RoMGDeDu7o5Vq1aVOp+6\nvgblya+ur0F+fj58fHzQuHFjuLm5Ydq0aaXOp67bvjz55dr2nIKKioq4unXrco8ePeIKCws5Dw8P\n7vbt2yXmOXjwINepUyeO4zju4sWLnI+Pj6KrVZry5D99+jQXEBDAKGHZ/vnnHy42NpZzd3cvdbo6\nb3uO+3x+dd72HMdxz5494+Li4jiO47js7GzOxcVFo77/y5NfnV+D3NxcjuM4TiKRcD4+Plx0dHSJ\n6eq87Tnu8/nl2fYKXyYgJiYGzs7OcHJygo6ODvr06YPw8PAS80RERGDQoEEAAB8fH2RkZOD58+eK\nrlopypMfUN8jeVq0aAEzM7NPTlfnbQ98Pj+gvtseAKytrdG4cWMAgLGxMVxdXfH06dMS86jza1Ce\n/ID6vgaGhoYAgMLCQkilUpibm5eYrs7bHvh8fqDi217hUk9NTS1xfLqdnR1SU1M/O8+TJ08UXbVS\nlCe/QCDA+fPn4eHhgc6dO+P27duVHVNu6rzty0OTtn1SUhLi4uLg4+NT4nFNeQ0+lV+dXwOZTIbG\njRvDysoKrVu3hpubW4np6r7tP5dfnm0v9xul/11peXz420ZdTkwqTw5PT0+kpKTA0NAQhw8fRlBQ\nEBISEiohnXKo67YvD03Z9jk5OQgODsbKlSthbGz80XR1fw3Kyq/Or4GWlhauXbuGzMxMdOjQAVFR\nURCLxSXmUedt/7n88mx7hffUbW1tkZKSUvx1SkoK7OzsypznyZMnanNZgfLkNzExKf4zqVOnTpBI\nJEhPT6/UnPJS521fHpqw7SUSCXr27In+/fsjKCjoo+nq/hp8Lr8mvAbVqlVDly5dcOXKlRKPq/u2\nf+9T+eXZ9gqXupeXF+7fv4+kpCQUFhZi586dCAwMLDFPYGAgwsLCAAAXL15E9erVYWVlpeiqlaI8\n+Z8/f1782z4mJgYcx5U69qWO1Hnbl4e6b3uO4zBs2DC4ublh/Pjxpc6jzq9BefKr62vw6tUrZGRk\nAADy8vJw/PhxiESiEvOo87YvT355tr3Cwy9CoRBr1qxBhw4dIJVKMWzYMLi6umL9+vUAgP/973/o\n3LkzDh06BGdnZxgZGWHz5s2KrlZpypN/z549WLduHYRCIQwNDbFjxw7Gqf/Vt29fnDlzBq9evYK9\nvT1mz54NiUQCQP23PfD5/Oq87QHg3Llz2Lp1Kxo1alT8A7lgwQIkJycDUP/XoDz51fU1ePbsGQYN\nGgSZTAaZTIYBAwagbdu2GtM95ckvz7aX+zIBhBBC1A/d+YgQQniESp0QQniESp0QQniESp0QQniE\nSp0QQniESp0QQnjk/wB0Qz82C0cRZwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10668b390>"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 7"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "t = np.random.rand(10000)\n",
      "moms = pd.DataFrame([[.5, t.mean()], [1 / np.sqrt(12), t.std(ddof=1)]],\n",
      "                    columns=['Theory', 'Data'],\n",
      "                    index=['Mean', 'Std'])\n",
      "moms"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>Theory</th>\n",
        "      <th>Data</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>Mean</th>\n",
        "      <td> 0.500000</td>\n",
        "      <td> 0.499425</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>Std</th>\n",
        "      <td> 0.288675</td>\n",
        "      <td> 0.287874</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "        Theory      Data\n",
        "Mean  0.500000  0.499425\n",
        "Std   0.288675  0.287874"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.random.rand(50).mean()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "0.51812069560161922"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 8"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.linalg import inv, lstsq, norm\n",
      "A = np.random.randn(500, 5)\n",
      "b = np.random.randn(500)\n",
      "eq = np.dot(inv(A.T.dot(A)).dot(A.T), b)\n",
      "least = lstsq(A, b)[0]\n",
      "diff = norm(eq - least)\n",
      "print('Error is %.3e' % diff)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Error is 1.724e-16\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Lab 5"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 1"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "# First matrix\n",
      "def flatDiags(n):\n",
      "    \"\"\"\n",
      "    Create n x n upper triangular matrix with values\n",
      "    on off diagonals equal to diag offset\n",
      "    \"\"\"\n",
      "    x = np.zeros((n, n))\n",
      "    for i in range(n):\n",
      "        x +=  np.diagflat(np.ones(n - i) * (i + 1), i)\n",
      "    return x\n",
      "\n",
      "flatDiags(5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "array([[ 1.,  2.,  3.,  4.,  5.],\n",
        "       [ 0.,  1.,  2.,  3.,  4.],\n",
        "       [ 0.,  0.,  1.,  2.,  3.],\n",
        "       [ 0.,  0.,  0.,  1.,  2.],\n",
        "       [ 0.,  0.,  0.,  0.,  1.]])"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Second type\n",
      "def totalDiags(n):\n",
      "    x = np.zeros((n,n))\n",
      "    for i in range(n):\n",
      "        x+= np.diagflat(np.ones(n-i)/float((i+1)),i) + np.diagflat(np.ones(n-i)/float((i+1)),-i)\n",
      "    x -= np.diag(np.ones(n),0)\n",
      "    return x\n",
      "\n",
      "totalDiags(5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "array([[ 1.        ,  0.5       ,  0.33333333,  0.25      ,  0.2       ],\n",
        "       [ 0.5       ,  1.        ,  0.5       ,  0.33333333,  0.25      ],\n",
        "       [ 0.33333333,  0.5       ,  1.        ,  0.5       ,  0.33333333],\n",
        "       [ 0.25      ,  0.33333333,  0.5       ,  1.        ,  0.5       ],\n",
        "       [ 0.2       ,  0.25      ,  0.33333333,  0.5       ,  1.        ]])"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 2"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.linalg as spla\n",
      "\n",
      "# First kind with toeplitz and triu\n",
      "def toep_triu(n):\n",
      "    return spla.triu(spla.toeplitz(np.arange(1, n + 1, 1)))\n",
      "\n",
      "toep_triu(5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "array([[1, 2, 3, 4, 5],\n",
        "       [0, 1, 2, 3, 4],\n",
        "       [0, 0, 1, 2, 3],\n",
        "       [0, 0, 0, 1, 2],\n",
        "       [0, 0, 0, 0, 1]])"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Second kind with toeplitz\n",
      "def toep_fracs(n):\n",
      "    return spla.toeplitz(1. / np.arange(1, n + 1))\n",
      "\n",
      "toep_fracs(5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 23,
       "text": [
        "array([[ 1.        ,  0.5       ,  0.33333333,  0.25      ,  0.2       ],\n",
        "       [ 0.5       ,  1.        ,  0.5       ,  0.33333333,  0.25      ],\n",
        "       [ 0.33333333,  0.5       ,  1.        ,  0.5       ,  0.33333333],\n",
        "       [ 0.25      ,  0.33333333,  0.5       ,  1.        ,  0.5       ],\n",
        "       [ 0.2       ,  0.25      ,  0.33333333,  0.5       ,  1.        ]])"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def ascending_diag(n):\n",
      "    return np.diagflat(np.arange(1, n+1, 1))\n",
      "\n",
      "ascending_diag(5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 24,
       "text": [
        "array([[1, 0, 0, 0, 0],\n",
        "       [0, 2, 0, 0, 0],\n",
        "       [0, 0, 3, 0, 0],\n",
        "       [0, 0, 0, 4, 0],\n",
        "       [0, 0, 0, 0, 5]])"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 3"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def second_der(n):\n",
      "    r1 = np.zeros(n)\n",
      "    r1[0] = 2\n",
      "    r1[1] = -1\n",
      "    return spla.toeplitz(r1)\n",
      "\n",
      "second_der(5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 25,
       "text": [
        "array([[ 2., -1.,  0.,  0.,  0.],\n",
        "       [-1.,  2., -1.,  0.,  0.],\n",
        "       [ 0., -1.,  2., -1.,  0.],\n",
        "       [ 0.,  0., -1.,  2., -1.],\n",
        "       [ 0.,  0.,  0., -1.,  2.]])"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 4"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.sparse as spar\n",
      "def second_der_sparse(n, format=None):\n",
      "    r1 = np.ones(n) * -1\n",
      "    r2 = np.ones(n) * 2\n",
      "    r3 = np.ones(n) * -1\n",
      "    \n",
      "    return spar.spdiags([r1, r2, r3], [-1, 0, 1], n, n, format=format)\n",
      "\n",
      "sparse_2nd = second_der_sparse(5)\n",
      "sparse_2nd"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 26,
       "text": [
        "<5x5 sparse matrix of type '<type 'numpy.float64'>'\n",
        "\twith 13 stored elements (3 diagonals) in DIAgonal format>"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sparse_2nd.todense()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 27,
       "text": [
        "matrix([[ 2., -1.,  0.,  0.,  0.],\n",
        "        [-1.,  2., -1.,  0.,  0.],\n",
        "        [ 0., -1.,  2., -1.,  0.],\n",
        "        [ 0.,  0., -1.,  2., -1.],\n",
        "        [ 0.,  0.,  0., -1.,  2.]])"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 5"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.sparse import linalg as sparla\n",
      "n_list = [50, 100, 300, 500, 1000, 3000]\n",
      "spar_times = []\n",
      "den_times = []\n",
      "for n in n_list:\n",
      "    A_spar = second_der_sparse(n)\n",
      "    A = second_der(n)\n",
      "    b = np.random.randn(n)\n",
      "    \n",
      "    # time sparse\n",
      "    t1 = time()\n",
      "    sparla.spsolve(A_spar, b)\n",
      "    spar_times.append(time() - t1)\n",
      "    \n",
      "    # time sparse\n",
      "    t1 = time()\n",
      "    spla.solve(A, b)\n",
      "    den_times.append(time() - t1)\n",
      "    \n",
      "times = pd.DataFrame([spar_times, den_times], \n",
      "                      columns=n_list, index=['Sparse', 'Dense']).T\n",
      "times"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/Users/spencerlyon2/anaconda/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:59: SparseEfficiencyWarning: spsolve requires CSC or CSR matrix format\n",
        "  warn('spsolve requires CSC or CSR matrix format', SparseEfficiencyWarning)\n"
       ]
      },
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>Sparse</th>\n",
        "      <th>Dense</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>50  </th>\n",
        "      <td> 0.039803</td>\n",
        "      <td> 0.000579</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>100 </th>\n",
        "      <td> 0.001207</td>\n",
        "      <td> 0.000641</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>300 </th>\n",
        "      <td> 0.001133</td>\n",
        "      <td> 0.003571</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>500 </th>\n",
        "      <td> 0.001297</td>\n",
        "      <td> 0.008589</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>1000</th>\n",
        "      <td> 0.001688</td>\n",
        "      <td> 0.074200</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>3000</th>\n",
        "      <td> 0.004103</td>\n",
        "      <td> 1.710525</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "output_type": "pyout",
       "prompt_number": 28,
       "text": [
        "        Sparse     Dense\n",
        "50    0.039803  0.000579\n",
        "100   0.001207  0.000641\n",
        "300   0.001133  0.003571\n",
        "500   0.001297  0.008589\n",
        "1000  0.001688  0.074200\n",
        "3000  0.004103  1.710525"
       ]
      }
     ],
     "prompt_number": 28
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 6"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vals, _ = spla.eig(second_der_sparse(1200).todense())\n",
      "a = vals.min()\n",
      "a * 1200 ** 2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 29,
       "text": [
        "(9.8531699800200343+0j)"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# This looks like pi ^ 2\n",
      "np.pi ** 2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 30,
       "text": [
        "9.869604401089358"
       ]
      }
     ],
     "prompt_number": 30
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Lab 6"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 1"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import numpy.linalg as npla\n",
      "\n",
      "A = np.array([[.75, .5], [.25, .5]])\n",
      "x = np.array([1, 0])\n",
      "ans = npla.matrix_power(A, 3).dot(x)\n",
      "\n",
      "print('The probablility of making the 3rd given he '+\n",
      "      'made the first is %.3f%%\\n' % ans[0])\n",
      "\n",
      "long_run = npla.matrix_power(A, 100)\n",
      "print('Long run matrix:\\n %s' % long_run)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The probablility of making the 3rd given he made the first is 0.672%\n",
        "\n",
        "Long run matrix:\n",
        " [[ 0.66666667  0.66666667]\n",
        " [ 0.33333333  0.33333333]]\n"
       ]
      }
     ],
     "prompt_number": 31
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can see that the long run probability of making a free throw is 0.667%"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 2"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# The matrix implied by the transition diagram\n",
      "A2 = np.array([[.25, 1 / 3., .5],\n",
      "               [.25, 1 / 3., 1 / 3.],\n",
      "               [.5, 1 / 3., 1 / 6.]])\n",
      "\n",
      "x = np.array([1, 0, 0])\n",
      "ans = npla.matrix_power(A2, 2).dot(x)\n",
      "print('The probability of being in state B two periods after '+\n",
      "      'being in state A is %.3f%%\\n' % ans[1])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The probability of being in state B two periods after being in state A is 0.312%\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "npla.matrix_power(A2, 50)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 33,
       "text": [
        "array([[ 0.35955056,  0.35955056,  0.35955056],\n",
        "       [ 0.30337079,  0.30337079,  0.30337079],\n",
        "       [ 0.33707865,  0.33707865,  0.33707865]])"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Because all the columns of the above matrix are the same, we say that the stable fixed piont is [0.3596, 0.3034, 0.3371]"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 3"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "A3 = np.array([[0, 0, 1, 0, 1, 0, 1],\n",
      "               [1, 0, 0, 0, 0, 1, 0],\n",
      "               [0, 0, 0, 0, 0, 1 ,0],\n",
      "               [1, 0, 0, 0, 1, 0, 0],\n",
      "               [0, 0, 0, 1, 0, 0, 0],\n",
      "               [0, 0, 1, 0, 0, 0, 1],\n",
      "               [0, 1, 0, 0, 0, 0, 0]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "len_5 = npla.matrix_power(A3, 5)\n",
      "len_7 = npla.matrix_power(A3, 7)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 35
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "len_5"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 36,
       "text": [
        "array([[2, 3, 1, 2, 1, 4, 1],\n",
        "       [0, 2, 5, 1, 3, 2, 5],\n",
        "       [0, 1, 2, 0, 1, 1, 2],\n",
        "       [1, 2, 3, 2, 3, 2, 3],\n",
        "       [2, 0, 2, 1, 2, 1, 2],\n",
        "       [1, 2, 1, 1, 0, 3, 1],\n",
        "       [3, 0, 2, 0, 1, 2, 2]])"
       ]
      }
     ],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.where(len_5 == len_5.max())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 37,
       "text": [
        "(array([1, 1]), array([2, 6]))"
       ]
      }
     ],
     "prompt_number": 37
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There max number of 5 node paths between nodes is 5. They are from node $2 \\rightarrow 3$ and from $2 \\rightarrow 7$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "len_7"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 38,
       "text": [
        "array([[ 2,  6,  9,  4,  6,  7,  9],\n",
        "       [ 8,  2, 10,  1,  6,  7, 10],\n",
        "       [ 3,  1,  4,  0,  2,  3,  4],\n",
        "       [ 6,  3,  9,  3,  7,  6,  9],\n",
        "       [ 4,  3,  3,  3,  3,  5,  3],\n",
        "       [ 1,  4,  6,  2,  3,  5,  6],\n",
        "       [ 3,  5,  2,  3,  1,  7,  2]])"
       ]
      }
     ],
     "prompt_number": 38
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There are no paths from node 3 to node 4"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bucky = np.genfromtxt('../Data/bucky.csv', delimiter=',')\n",
      "print('number of connections in bucky: %i' % \n",
      "      np.count_nonzero(bucky))\n",
      "msg = 'number %i-node connections in bucky: %i'\n",
      "print(msg % (2, np.count_nonzero(npla.matrix_power(bucky, 2))))\n",
      "print(msg % (3, np.count_nonzero(npla.matrix_power(bucky, 3))))\n",
      "print(msg % (4, np.count_nonzero(npla.matrix_power(bucky, 4))))\n",
      "\n",
      "done = True\n",
      "i = 4  # Start where we left off above\n",
      "while done:\n",
      "    i += 1\n",
      "    non_zeros = np.count_nonzero(npla.matrix_power(bucky, i))\n",
      "    done = non_zeros != bucky.size\n",
      "\n",
      "print('At path length %i, all atoms are connected' % i)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "number of connections in bucky: 180\n",
        "number 2-node connections in bucky: 420\n",
        "number 3-node connections in bucky: 780\n",
        "number 4-node connections in bucky: 1380\n",
        "At path length 9, all atoms are connected\n"
       ]
      }
     ],
     "prompt_number": 39
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def spy_plot(n):\n",
      "    plt.figure()\n",
      "    plt.spy(npla.matrix_power(bucky, n))\n",
      "    plt.title('Path Length %i' % n)\n",
      "\n",
      "spy_plot(1)\n",
      "spy_plot(2)\n",
      "spy_plot(4)\n",
      "spy_plot(10)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAEHCAYAAABlS0A3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF1NJREFUeJzt3X9QFOf9B/D3aZgqDeqhuQU1ydlEJCDCWUdnbH6s0gPb\nCMHSUDMtw6jT6XTqpGamzaiZBpK2cramGdPmn3b645pMMTR/KBhDZKLHJLEpY4pjUpPSGghUORpz\nR+SHBIH9/pEvK4dwe8vd3u7d837N7Ax3t7f7uT0+9zyf3Wd3bYqiKCAiYcwxOwAiii8mPZFgmPRE\ngmHSEwmGSU8kGCY9kWCY9Elgzpw5+PDDD80OQxdZlvH73//e7DCExKQ3idPpRGpqKtLS0pCRkYEd\nO3ZgcHBQ833RJktNTQ0qKytn/f5YrdNms8Fms0X0/uvXr+Ob3/wmVqxYgTlz5qClpcWIMIXBpDeJ\nzWbD8ePH0d/fj3/84x84e/Ysfvazn0X0vmjXm4juv/9+vPjii8jIyEjYz2AVTHoLWLp0KbZs2YL3\n3nsPfX192Lp1KxwOB9LT01FSUoJLly4BAJ544gm88cYb2L17N9LS0vDoo4+qy2hubkZWVhbsdjt2\n794947rCDcB8++23sXHjRtjtdhQUFIS0qLIs48knn8S9996LBQsWoLi4GJ988on6+p///Gfceeed\nWLJkCX7605/C6XTi9ddfR1NTE2pra/HSSy8hLS0NLpdLfU9nZ+eMy5ssJSUFjz76KL7yla9g7ty5\n2huUwmLSm2giAbu7u/Hqq69i7dq1GB8fx65du9DV1YWuri7Mnz9fTeKf//znuO+++/D888+jv78f\nzz33nLqsV155BWfPnsX58+dRX1+P1157TVcsly5dwtatW/Hkk08iGAzi0KFDKC8vD0nEuro6/OlP\nf8L//vc/jIyM4NChQwCACxcu4Ac/+AHq6urQ09ODq1ev4vLly7DZbNiyZQv279+P7du3o7+/H21t\nbepn/8tf/jLt8shYTHqTKIqCsrIy2O123HfffZBlGfv370d6ejq2bduGefPm4dZbb8X+/ftvqmGn\na6337t2LBQsW4Pbbb8emTZtw7tw5XfG8+OKL+PrXv44tW7YAAL761a9i3bp1eOWVVwB8Xhbs2LED\nd999N+bNm4eKigp1HS+//DJKS0uxceNGpKSk4Omnnw7pgiuKclPMNpsNO3funHZ5ZKxbzA5AVDab\nDceOHcPmzZtDnh8aGsJjjz2G1157DcFgEAAwMDAARVHURJqups3IyFD/Tk1NxcDAgK54PvroI/z1\nr39FY2Oj+tzo6GhIfJPXMX/+fHUdly9fxvLly0NeW7x4seY6Z1oeGYtJbzHPPPMM2tvb0draCofD\ngXPnzmHt2rVq0hu1I++OO+5AZWUlfvvb3+pe5tKlS/Gvf/1LfXzt2rWQsoA73qyF3XuLGRgYwPz5\n87Fw4UIEAgE89dRTIa9LkoSLFy+GXUa4nXWKomB8fByfffYZhoeHMTw8jM8++wzf+c530NjYiJMn\nT2JsbAzDw8Pw+XzqTsRwyy0vL0djYyP+9re/YWRkBDU1NSHzZmRkoLOz86b36zmreyLeqX+Tfkx6\ni9mzZw+uXbuGJUuWYOPGjfja174W0lL+8Ic/xMsvv4z09HTs2bNn2mWE6xHYbDbU1dVh/vz5SE1N\nRWpqKlauXInly5fj2LFjOHDgABwOB+644w4888wzIYk5eZmT15Gbm4tf//rX2L59O5YuXYq0tDQ4\nHA584QtfAAA8/PDDAIDFixdj3bp1msubzqpVq5CamorLly+juLgYX/ziF9HV1TXj/DQzGy+iQbE2\nMDAAu92O//znP7jzzjvNDoemYEtPMdHY2IihoSEMDg7iRz/6EdasWcOEtygmPcVEQ0MDli1bhmXL\nluHixYs4cuSI2SHRDNi9JxKMoS19U1MTsrOzsXLlShw8eNDIVc1o586dkCQJeXl56nOBQAButxtZ\nWVkoKipCX19fXGPq7u7Gpk2bkJubi9WrV6sj68yMa3h4GBs2bEBBQQFycnKwb98+02OaMDY2BpfL\nhZKSEsvE5HQ6sWbNGrhcLqxfv94ycUXCsKQfGxvD7t270dTUhAsXLqCurg7vv/++Uaub0Y4dO9DU\n1BTynMfjgdvtRnt7OwoLC+HxeOIaU0pKCp599ln885//xNtvv43nn38e77//vqlxzZs3D6dPn8a5\nc+dw/vx5nD59Gm+++abp2woADh8+jJycHHXvvhVistls8Pl8aGtrQ2trq2XiiohikDNnzijFxcXq\n49raWqW2ttao1YXV0dGhrF69Wn28atUqxe/3K4qiKD09PcqqVatMiWvCQw89pDQ3N1smrsHBQWXd\nunXKe++9Z3pM3d3dSmFhoXLq1Cll69atiqJY4/tzOp3KlStXQp6zQlyRMKylv3TpEm6//Xb18fLl\ny0MGepipt7cXkiQB+HywS29vr2mxdHZ2oq2tDRs2bDA9rvHxcRQUFECSJLX8MDumxx57DL/85S8x\nZ86Nf1WzYwI+b+knzk/43e9+Z5m4ImHYMNxEGXoZi6GtszUwMIDy8nIcPnwYaWlppsc1Z84cnDt3\nDp9++imKi4tx+vRpU2M6fvw4HA4HXC4XfD7ftPOY9f299dZbyMzMxMcffwy3243s7GxLxBUJw1r6\nZcuWobu7W33c3d0dclKGmSRJgt/vBwD09PTA4XDEPYbr16+jvLwclZWVKCsrs0xcALBw4UI8+OCD\neOedd0yN6cyZM2hoaMCKFSvwyCOP4NSpU6isrLTEdsrMzAQA3Hbbbdi2bRtaW1stEVckDEv6devW\n4d///jc6OzsxMjKCl156CaWlpUatTpfS0lJ4vV4AgNfrVZMuXhRFwa5du5CTkxMylNbMuK5cuaLu\nbb527Rqam5vhcrlMjenAgQPo7u5GR0cHjhw5gs2bN+OFF14w/fsbGhpCf38/AGBwcBAnT55EXl6e\n6XFFzMgdBidOnFCysrKUu+66Szlw4ICRq5rR9u3blczMTCUlJUVZvny58oc//EH55JNPlMLCQmXl\nypWK2+1WgsFgXGN64403FJvNpuTn5ysFBQVKQUGB8uqrr5oa1/nz5xWXy6Xk5+creXl5yi9+8QtF\nURTTt9UEn8+nlJSUWCKmDz/8UMnPz1fy8/OV3Nxc9X/b7LgixcE5RILhMFwiwcw66a0w2o6IZmE2\nNcHo6Khy1113KR0dHcrIyIiSn5+vXLhwIWQeAJw4cTJpCmdWLX1rayvuvvtuOJ1OpKSkYPv27Th2\n7NhN81VXV+OBBx5AdXX1Tcd8lf+/WOJMUzha79Waqquro15GrCfGlLgxmR3X6dOnUV1drU5aZjU4\nZ7rRdn//+99vmq+mpkadiMgYsixDlmX18dRLrE01q5beqiONiEjbrFr6SEfbTfw4TPzyKIpy02sT\nJr823eNYCReTmSb/UlsFY4qcVeOazqyO04+OjmLVqlV4/fXXsXTpUqxfvx51dXW45557bix4mt6A\nnqQ3ilnrJYoXm80W9v96Vi39Lbfcgt/85jcoLi7G2NgYdu3aFZLwRGRdho3Im+7XZuqtjmZ6bbrX\nrSgRY6bkp9XSc0QekWCY9ESCYdITCSauN7DUs/c+EerlWMYUbn8HUSyxpScSDJOeSDBMeiLBGFrT\nh6vLE7GGN1KsPq/o25G0saUnEgyTnkgwTHoiwcT1OH04empPvXWrSMfAk/3zUfTY0hMJhklPJBhD\nu/fxGqZq1lV4kg0P94mBLT2RYJj0RIJh0hMJxjKH7LQk22m5VpSI24nftX5s6YkEw6QnEgyTnkgw\nCVPTT6a3hhdpGG4iiqYu5/epH1t6IsEw6YkEw6QnEkzC1PTh6vJojtPzOK/5uM3jiy09kWDCJv3O\nnTshSRLy8vLU5wKBANxuN7KyslBUVIS+vj7DgySi2Amb9Dt27EBTU1PIcx6PB263G+3t7SgsLITH\n44l4ZTabTZ30UhRFnfTMG+vTe6P5DERWoHmr6s7OTpSUlODdd98FAGRnZ6OlpQWSJMHv90OWZXzw\nwQc3L1jnrarNwvqfko3Wrap178jr7e2FJEkAAEmS0NvbO+O8NTU16t+yLOtdFRFFwOfzwefzRTy/\n7pbebrcjGAyqr6enpyMQCNy8YLb0RKbQaul1772f6NYDQE9PDxwOR8TvNaLOjsTUWnzyNLX+D1e3\nG7WvgCiedCd9aWkpvF4vAMDr9aKsrCzmQRGRccJ27x955BG0tLTgypUrkCQJTz/9NB566CFUVFSg\nq6sLTqcT9fX1WLRo0c0L1uhixFO4Pe28AAclG63c06zpjVpxPDHpSSQx33ufiKLZOWfFnY+Jjj+s\n5uIwXCLBMOmJBMOkJxKMEDV9OHp25LEWvYGXuEpcbOmJBMOkJxKMcN17vXe41XPITqTDe8n++ZIZ\nW3oiwTDpiQTDpCcSjGk1vVmHv/Sux6i4ePiPzMKWnkgwTHoiwTDpiQRjWk2fiDWs1mWv9QzZTcTP\nT8mBLT2RYJj0RIIRbhhuNKIZhpsMh+hEGmaczNjSEwmGSU8kGCY9kWBY08eQnkN2iVjjJ0KMVmS1\n75otPZFgmPREgmHSEwmGNX0UwtVq0dTwVqsBKTpW+/7Y0hMJJmzSd3d3Y9OmTcjNzcXq1avx3HPP\nAQACgQDcbjeysrJQVFSEvr6+uARLRNELe9dav98Pv9+PgoICDAwM4Mtf/jKOHj2KP/7xj1iyZAke\nf/xxHDx4EMFgEB6PJ3TBFrprrVGM6rKze0/R0Mq9sC19RkYGCgoKAAC33nor7rnnHly6dAkNDQ2o\nqqoCAFRVVeHo0aMxCXTylAgURQmZYjWvXom23chcEe/I6+zsRFtbGzZs2IDe3l5IkgQAkCQJvb29\n076npqZG/VuWZciyHFWwRHQzn88Hn88X8fxhu/cTBgYG8MADD+AnP/kJysrKYLfbEQwG1dfT09MR\nCARCF6yze88u7Q16twXPfqPJoureA8D169dRXl6OyspKlJWVAfi8dff7/QCAnp4eOByOqAM1svtr\nRVPLmcnT1G2hVfqItN0oemGTXlEU7Nq1Czk5OdizZ4/6fGlpKbxeLwDA6/WqPwZEZH1hu/dvvvkm\n7r//fqxZs0ZtYWpra7F+/XpUVFSgq6sLTqcT9fX1WLRoUeiCBdh7H41wO92S4eQcMo9W7kVU0xux\nYtEx6ckoWrnHYbgmieY4PX8ErM3q3w+H4RIJhklPJBgmPZFgWNNbkN4anoNzYi+autzq3wFbeiLB\nMOmJBMOkJxIMa3qLCFeX89Jb8ZfM24ktPZFgmPREgrFs9160bqkVDgmJts1FxZaeSDBMeiLBMOmJ\nBGPZmp71ZORidW4+t7kY2NITCYZJTyQYJj2RYCxb01PkYnXpLdFrelG2BVt6IsEw6YkEw6QnEgxr\n+iSn5zh9Moy9T+bLXMUKW3oiwTDpiQSTkN37ZOiGGkVr2yT7Ibtk+AxGY0tPJJiwST88PIwNGzag\noKAAOTk52LdvHwAgEAjA7XYjKysLRUVF6Ovri0uwRBQ9zbvWDg0NITU1FaOjo7j33ntx6NAhNDQ0\nYMmSJXj88cdx8OBBBINBeDye0AUbeNdadu9npnfbJFv3nrRzT7N7n5qaCgAYGRnB2NgY7HY7Ghoa\nUFVVBQCoqqrC0aNHYxRuZBRFCZnoBr3bRs+8NptNnfTMG8n8IjF722juyBsfH8fatWtx8eJFfP/7\n30dubi56e3shSRIAQJIk9Pb2Tvvempoa9W9ZliHLckyCJqIbfD4ffD5fxPNrdu8nfPrppyguLkZt\nbS2+8Y1vIBgMqq+lp6cjEAiELtjA7j2ZQ08pwBJsZkZvm6i79xMWLlyIBx98EO+88w4kSYLf7wcA\n9PT0wOFwRB8pEcVF2KS/cuWKumf+2rVraG5uhsvlQmlpKbxeLwDA6/WirKzM+EgTjNl1WyxM/QyT\n63+tz8f9LjMze9uE7d6/++67qKqqwvj4OMbHx1FZWYkf//jHCAQCqKioQFdXF5xOJ+rr67Fo0aLQ\nBQvevU+G7m2yj9NPVlq5F3FNH+sVJ7tkSAomfWLSyr2EHIYbT7M9jp0MSaBnyC5/BBIHh+ESCYZJ\nTyQYJj2RYFjTa2BteoOeHXl67rpD8cWWnkgwTHoiwTDpiQTDmp5mpOcOtzxOHxvx2I5s6YkEw6Qn\nEgyTnkgwrOlpRla5xqFI1/GLx+djS08kGCY9kWDYvZ+Ch55mJ5ouut678lB02NITCYZJTyQYJj2R\nYJKypo+mLmf9aIxoTsvldxJbbOmJBMOkJxIMk55IMElZ07MGtDa9NTwvvRVbbOmJBMOkJxIMk55I\nMElZ05P1hBt7z+P0MzNiW7ClJxJMREk/NjYGl8uFkpISAEAgEIDb7UZWVhaKiorUe9gTkfVFlPSH\nDx9GTk6O2tXweDxwu91ob29HYWEhPB6PoUFS4lMURZ1iyWazhUyxmtcqJm+3WG07zZr+v//9L06c\nOIEnnngCv/rVrwAADQ0NaGlpAQBUVVVBluVpE7+mpkb9W5ZlyLIck6CJ6Aafzwefzxfx/DZF4+fj\n4Ycfxv79+3H16lUcOnQIjY2NsNvtCAaDAD7/JUpPT1cfqwu22YTeAUOzp2fnlVHzJjKt3AvbvT9+\n/DgcDgdcLteMC0mkrhIRaXTvz5w5g4aGBpw4cQLDw8O4evUqKisrIUkS/H4/MjIy0NPTA4fDEa94\npyXKL7goommt9dyVR1RhW/oDBw6gu7sbHR0dOHLkCDZv3owXXngBpaWl8Hq9AACv14uysrK4BEtE\n0dN1nH7iV3Tv3r1obm5GVlYWTp06hb179xoSHBHFnuaOvFkvOI478ti9FwevuqNNK/eSYhiuiF9s\nrCTa3WOiOS03ET5fPHAYLpFgmPREgmHSEwkmKWp6mr1Er3P11Phag8gSfVtEii09kWCY9ESCYfee\nEoreO9wm2yG7WIxDYEtPJBgmPZFgmPREgmFNTwnFyLpcT/1v1hh/Xg2XiHRj0hMJhklPJBjW9IIR\n7Xzz2V56S+94gETClp5IMEx6IsGwe5+AoumiJ3K3NNb0DNlNprKILT2RYJj0RIJh0hMJhjV9Akrk\netLK9ByyC3cVHqt/P2zpiQTDpCcSDJOeSDCs6UlYeu5wa9Xj9LO5HBhbeiLBaLb0TqcTCxYswNy5\nc5GSkoLW1lYEAgF861vfwkcffQSn04n6+nosWrQoHvESUZQ0W3qbzQafz4e2tja0trYCADweD9xu\nN9rb21FYWAiPx2N4oEQUG5q3ql6xYgXOnj2LxYsXq89lZ2ejpaUFkiTB7/dDlmV88MEHoQu22VBd\nXa0+lmUZsizHNnqiONFT0+ut/6O9TLfP54PP51MfP/XUU+Hj00r6L33pS1i4cCHmzp2L733ve/ju\nd78Lu92OYDCoBpmenq4+nvxBrLKzgyhaVk766ZYX1f3p33rrLWRmZuLjjz+G2+1Gdnb2TSvQukcY\nEVmHZk2fmZkJALjtttuwbds2tLa2qt16AOjp6YHD4TA2SsFN/LDyB9Y8iqKETFNN/n6mzqv1/YVb\nrhHCJv3Q0BD6+/sBAIODgzh58iTy8vJQWloKr9cLAPB6vSgrKzM+UiKKibA1fUdHB7Zt2wYAGB0d\nxbe//W3s27cPgUAAFRUV6OrqmvGQHWv62LHqwBC6wUoX4NDKPc0deUatmCLHpLe+REp6DsONE17i\nKrlFc1puvL9fDsMlEgyTnkgwTHoiwbCmj5NErMvNrj0Tld4aPtYj8rSwpScSDJOeSDBMeiLBsKan\nGbGGj1y4ujxex+kjPS+DLT2RYJj0RIJh954MIdrhPqOGVRuxHdnSEwmGSU8kGCY9kWBY05MhYlnD\nx3uYqpmiOS134m+tQ3ds6YkEw6QnEgyTnkgwrOnJ8owapmrF/QPxOC2XLT2RYJj0RIJh0hMJhjU9\nCcOosQOxXna45YZbL0+tJaJpGZ70k++bbRWMKTKMKXJWjWs6THqLYEyRsUpMU+9MG01cU+9qG+kd\nbqfeETfSO9+ye08kGCY9kWAMvWstEZnDlLvWWnGIIxGxe08kHCY9kWCY9ESCYdITCYZJTyQYJj2R\nYP4PltjHk8Sjpz8AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106660d50>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAEHCAYAAABlS0A3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGVVJREFUeJzt3X9QFOcZB/DvmTAVRlTQ3IGa5NpENCDCWUdnbNKe0gOb\nCMEyoWZShkGnvyZOa2Y6GTVTQ5NWzsa0o23+aZu213SKoflDwRgiEz2aaC2TFMempqUx/CoBGnNH\n5IeEANs/Uk7uPG5vb3dv9+79fmZuBri93eeWe+59n91337VIkiSBiIQxz+gAiCi+mPREgmHSEwmG\nSU8kGCY9kWCY9ESCYdIngXnz5uG9994zOgxFnE4nnn/+eaPDEBKT3iB2ux1paWlIT09HVlYWampq\nMDo6Kvs6tclSW1uLqqqqmF+v1TYtFgssFktUr79w4QJcLheWLFkCq9WKyspKDAwM6BGqEJj0BrFY\nLDh58iSGh4fxt7/9DW+++SZ+9KMfRfU6tdtNNENDQ/j2t7+N7u5udHd3Iz09HTU1NUaHlbCY9Caw\nbNkybN26FW+//TaGhoawbds2WK1WZGZmorS0FH19fQCAJ554Aq+//jp2796N9PR0fPe73w2so6Wl\nBTk5OcjIyMDu3bvn3FakAZgXLlzApk2bkJGRgcLCQrS2tgaeczqdOHDgAO69914sXLgQJSUl+PDD\nDwPP//73v8edd96JpUuX4umnn4bdbsdrr72G5uZm1NXV4cUXX0R6ejocDkfgNV1dXXOub7atW7ei\noqICCxYsQGpqKh599FGcO3dOfsdSWEx6A80kYG9vL1555RWsW7cO09PT2LVrF3p6etDT04PU1NRA\nEv/4xz/Gfffdh+eeew7Dw8M4evRoYF0vv/wy3nzzTVy6dAkNDQ149dVXFcXS19eHbdu24cCBA/D7\n/Th8+DAqKiqCErG+vh6/+93v8N///hcTExM4fPgwAODy5ct49NFHUV9fj/7+fly7dg3vv/8+LBYL\ntm7div3792PHjh0YHh5Ge3t74L3/8Y9/DLs+OX/+85+xZs0aRe+PbmDSG0SSJJSXlyMjIwP33Xcf\nnE4n9u/fj8zMTGzfvh3z58/HggULsH///qAWd+a1ofbu3YuFCxfi9ttvx+bNm3Hx4kVF8fzhD3/A\n/fffj61btwIAvvzlL2P9+vV4+eWXAXxaFtTU1ODuu+/G/PnzUVlZGdjGSy+9hLKyMmzatAkpKSl4\n6qmngsoISZJuitlisWDnzp1h1xfJpUuX8PTTT+OZZ55R9P7ohluNDkBUFosFJ06cwJYtW4L+PjY2\nhsceewyvvvoq/H4/AGBkZASSJAUSKVxdnpWVFfg5LS0NIyMjiuLp7u7Gn/70JzQ1NQX+Njk5GRTf\n7G2kpqYGtvH+++9jxYoVQc8tWbJEdptzrW8u7777Lu6//34cPXoUX/jCF+TfFIXFpDeZZ599Fh0d\nHWhra4PVasXFixexbt26QNLrdSDvjjvuQFVVFX75y18qXueyZcvwr3/9K/D79evXg8oCLQ4ednd3\nw+Vy4cCBA3jkkUdUr09k7N6bzMjICFJTU7Fo0SL4fD788Ic/DHreZrPhypUrEdcR6WCdJEmYnp7G\nxx9/jPHxcYyPj+Pjjz/G17/+dTQ1NeH06dOYmprC+Pg4vF5v4CBipPVWVFSgqakJf/nLXzAxMYHa\n2tqgZbOystDV1XXT66O9qruvrw9btmzB7t278c1vfjOq19DcmPQms2fPHly/fh1Lly7Fpk2b8JWv\nfCWopfze976Hl156CZmZmdizZ0/YdUTqEVgsFtTX1yM1NRVpaWlIS0vDypUrsWLFCpw4cQIHDx6E\n1WrFHXfcgWeffTYoMWevc/Y28vLy8POf/xw7duzAsmXLkJ6eDqvVis985jMAgIceeggAsGTJEqxf\nv152faF+/etfo7OzE7W1tUhPT0d6ejoWLlwYcT/S3CycRIO0NjIygoyMDLz77ru48847jQ6HQrCl\nJ000NTVhbGwMo6Oj+P73v4+1a9cy4U2KSU+aaGxsxPLly7F8+XJcuXIFx44dMzokmgO790SC0bWl\nb25uxurVq7Fy5UocOnRIz03NaefOnbDZbMjPzw/8zefzweVyIScnB8XFxRgaGoprTL29vdi8eTPy\n8vKwZs2awMg6I+MaHx/Hxo0bUVhYiNzcXOzbt8/wmGZMTU3B4XCgtLTUNDHZ7XasXbsWDocDGzZs\nME1c0dAt6aemprB79240Nzfj8uXLqK+vxzvvvKPX5uZUU1OD5ubmoL+53W64XC50dHSgqKgIbrc7\nrjGlpKTgZz/7Gf7xj3/gwoULeO655/DOO+8YGtf8+fNx9uxZXLx4EZcuXcLZs2fxxhtvGL6vAODI\nkSPIzc0NHN03Q0wWiwVerxft7e1oa2szTVxRkXRy/vx5qaSkJPB7XV2dVFdXp9fmIurs7JTWrFkT\n+H3VqlXSwMCAJEmS1N/fL61atcqQuGY8+OCDUktLi2niGh0dldavXy+9/fbbhsfU29srFRUVSWfO\nnJG2bdsmSZI5/n92u126evVq0N/MEFc0dGvp+/r6cPvttwd+X7FiRdBADyMNDg7CZrMB+HSwy+Dg\noGGxdHV1ob29HRs3bjQ8runpaRQWFsJmswXKD6Njeuyxx/DMM89g3rwbH1WjYwI+belnrk/41a9+\nZZq4oqHbMNxEuW5bi6GtsRoZGUFFRQWOHDmC9PR0w+OaN28eLl68iI8++gglJSU4e/asoTGdPHkS\nVqsVDocDXq837DJG/f/OnTuH7OxsfPDBB3C5XFi9erUp4oqGbi398uXL0dvbG/i9t7c36KIMI9ls\ntsDMK/39/bBarXGP4ZNPPkFFRQWqqqpQXl5umrgAYNGiRXjggQfw1ltvGRrT+fPn0djYiM9+9rN4\n+OGHcebMGVRVVZliP2VnZwMAbrvtNmzfvh1tbW2miCsauiX9+vXr8e9//xtdXV2YmJjAiy++iLKy\nMr02p0hZWRk8Hg8AwOPxBJIuXiRJwq5du5Cbmxs0lNbIuK5evRo42nz9+nW0tLTA4XAYGtPBgwfR\n29uLzs5OHDt2DFu2bMELL7xg+P9vbGwMw8PDAIDR0VGcPn0a+fn5hscVNT0PGJw6dUrKycmR7rrr\nLungwYN6bmpOO3bskLKzs6WUlBRpxYoV0m9+8xvpww8/lIqKiqSVK1dKLpdL8vv9cY3p9ddflywW\ni1RQUCAVFhZKhYWF0iuvvGJoXJcuXZIcDodUUFAg5efnSz/5yU8kSZIM31czvF6vVFpaaoqY3nvv\nPamgoEAqKCiQ8vLyAp9to+OKFgfnEAmGw3CJBBNz0pthtB0RxSCWmmByclK66667pM7OTmliYkIq\nKCiQLl++HLQMAD744MOgRyQxtfRtbW24++67YbfbkZKSgh07duDEiROyr5P+P0FiNIcRZi8b+lCy\nrNy25F4br8eTTz6pyXq0fH9axWTG/ZRMcZ09exZPPvlk4CEnpsE54Ubb/fWvf41lVUSkktPphNPp\nDPweOsVaqJhaerOONCIieTElfbSj7Wa6HjPdkJmhiRaLRbZbOnvZ0Ifca8PFEdodCheT3JdZpJjU\nfhHO/qZWQ+m+iUdMWjJjTIB54wonpvP0k5OTWLVqFV577TUsW7YMGzZsQH19Pe65554bK/5/cgZt\nLOQGCHM9J0fthznSdiOtWy5GLeMiilW43Jstppr+1ltvxS9+8QuUlJRgamoKu3btCkp4IjIv3Ubk\nhWsVZ29KroVV0gKriUuvXoPceiP1GthjIDXkWnqOyCMSDJOeSDBMeiLBxPUGlkqO3utV1xp15F9J\nHErXy+MDpARbeiLBMOmJBMOkJxKMrjV9pLpd6Xn6SMsaxahz/EriUDPi0Cz7mbTFlp5IMEx6IsEw\n6YkEE9fz9LPpeb5cybbMUscmwrgESg5s6YkEw6QnEoyu3XujhqmGPh/pdFiydX/VDOFNxH1hlvIs\nkbClJxIMk55IMEx6IsHEdRiuGmouyxVpym41Q3iNqo+1HCrMGl8eW3oiwTDpiQTDpCcSjGHDcNXQ\nsq5jDXiDWWr4SM+zhlePLT2RYJj0RIJh0hMJxrDbWmm97kjbUTL1FmvE+OM+1xZva0VEQSIm/c6d\nO2Gz2ZCfnx/4m8/ng8vlQk5ODoqLizE0NKR7kESknYhJX1NTg+bm5qC/ud1uuFwudHR0oKioCG63\nO6YNWyyWoIcSkiRFfGhJSYyxvh/R6fn/o5vJ1vRdXV0oLS3F3//+dwDA6tWr0draCpvNhoGBATid\nTvzzn/+8ecUqb1WtJTU1faRlI22HH14yilxNr3hwzuDgIGw2GwDAZrNhcHAw9uiISDWv1wuv1xv1\n8opb+oyMDPj9/sDzmZmZ8Pl8N6+YLX2U0RFpS/OWfqZbn5WVhf7+flit1jmXjZRQSpNCyWv1ulRT\nblkmenR4is5Yik/ZlZWVwePxAAA8Hg/Ky8s1D4qI9BOxe//www+jtbUVV69ehc1mw1NPPYUHH3wQ\nlZWV6Onpgd1uR0NDAxYvXnzzisN0MRK9pQ/FFio2bOn1Jde913VEHpOewmHS60vzml4NNf9cNdNn\nazUMl+ffb9Byqm1+CcQXh+ESCYZJTyQYJj2RYBJyuiylItXlSg7kiT44R01dboaptulTbOmJBMOk\nJxJMXLv3Zuj+KrnDbejySmbsScYurFbvLxn2RSJjS08kGCY9kWCY9ESCMc0pO6NqYLPc1SXWZY1i\nxpgoOmzpiQTDpCcSDJOeSDCGXVprVN2q5SWhSpZVMh7ArPWyGcZZkHps6YkEw6QnEoxhp+zMMgxX\nyfJaDsNNhFl4lAxJZnd/bmbbV2zpiQTDpCcSDJOeSDCmGYabCJSchlMz46tZakAtpzAXidLjOfHe\nl2zpiQTDpCcSDJOeSDCG1fRG1zVaUHOHW63uuhNPifg/MoJZ/38z2NITCSZi0vf29mLz5s3Iy8vD\nmjVrcPToUQCAz+eDy+VCTk4OiouLMTQ0FJdgiUi9iHetHRgYwMDAAAoLCzEyMoLPf/7zOH78OH77\n299i6dKlePzxx3Ho0CH4/X643e7gFcvdLtdkXR4tqOmyK1k2GfaVSOL9/5PLvYgtfVZWFgoLCwEA\nCxYswD333IO+vj40NjaiuroaAFBdXY3jx49rGPKnLBZL0CMRSJI050NLSvZNIu7HZKPnZyEWUR/I\n6+rqQnt7OzZu3IjBwUHYbDYAgM1mw+DgYNjX1NbWBn52Op1wOp2qgiWim3m9Xni93qiXj9i9nzEy\nMoIvfelL+MEPfoDy8nJkZGTA7/cHns/MzITP5wtescrufbJ3abUqBcItH+m1ybYf6WZyuSfb0n/y\nySeoqKhAVVUVysvLAXzaug8MDCArKwv9/f2wWq2KA0uGJNdraK3SL8RIw2MTYT8mG7N/liPW9JIk\nYdeuXcjNzcWePXsCfy8rK4PH4wEAeDyewJcBEZlfxO79G2+8gS9+8YtYu3Zt4Nurrq4OGzZsQGVl\nJXp6emC329HQ0IDFixcHr1imi3FTICb/dgzHLC19tOul+DD6syxbWkdT0+ux4XDLz5YIH14mPYVj\n9GdZdU0fL4n4YVVzp1Ytz9PzVJz2tJw12egvgVAchkskGCY9kWCY9ESCMU1NLxotL8vl4BxtqKnL\nE2l6cLb0RIJh0hMJhklPJBjT1PRmr4P0FqkmTMapt8xIq7rc7PuULT2RYJj0RILRtXuvV1cy2buo\n8Xw/vGtNeMm8L9jSEwmGSU8kGCY9kWB0rem1vPR09vPJWG8pqa31ujY/2Y+VyBHl+AZbeiLBMOmJ\nBMOkJxKMaYbhhkrmmiocM0y9lQzTbml1vCOZj2+wpScSDJOeSDBMeiLBmLamjySZ661YaHVZbiLu\nRzXvL9n2RbTY0hMJhklPJBjTdu85G0xslNzhNnT5RNyPWn4WEuH9aoEtPZFgIib9+Pg4Nm7ciMLC\nQuTm5mLfvn0AAJ/PB5fLhZycHBQXF2NoaCguwRKRerJ3rR0bG0NaWhomJydx77334vDhw2hsbMTS\npUvx+OOP49ChQ/D7/XC73cErVnjX2psCY/deE3Kj7BK9ex8qGd6DWnK5J9u9T0tLAwBMTExgamoK\nGRkZaGxsRHV1NQCguroax48f1yjcGyRJmvMht6wZWCyWoIdRIu3HcKewZj8ivYfQ58z6fs1Abt/E\nez/KHsibnp7GunXrcOXKFXznO99BXl4eBgcHYbPZAAA2mw2Dg4NhX1tbWxv42el0wul0ahI0Ed3g\n9Xrh9XqjXl62ez/jo48+QklJCerq6vDVr34Vfr8/8FxmZiZ8Pl/wilV27xNdMnQz1Vyck4jvVy9a\n3rcg2u2p6t7PWLRoER544AG89dZbsNlsGBgYAAD09/fDarUqDoyIjBEx6a9evRo4Mn/9+nW0tLTA\n4XCgrKwMHo8HAODxeFBeXq5/pAZRUl/NXlZJfWxWkd6DXO2ciO9XL3KfBSXHr7QQsabv7+9HdXU1\npqenMT09jaqqKhQVFcHhcKCyshLPP/887HY7GhoadAmOiLQXdU2veMVJUtPHOmFlMp5WVPL+QiXi\n+9WL3p8Fudwz7TBco6gZpirSOW+lB6dEuYItGka/fw7DJRIMk55IMEx6IsGwpg+h1WwqRtdtetDq\nQKXS4x0iHQ+Ix7EgtvREgmHSEwmGSU8kGA7OoaipmeMg0rKi0ev4R7THPtjSEwmGSU8kGCY9kWB4\nnj5EMoyZ14ted7gV7WIdvd7PzHrl9idbeiLBMOmJBJOU3XstT4mINARUCT2H4fKy3Bs4DJeIVGPS\nEwmGSU8kmKSo6dXU5ck4l50R9Jw/cPbz/P+ox5aeSDBMeiLBMOmJBJMUNb2WdTlrRP0pPU+vx6Wn\ncsuahR4xs6UnEgyTnkgwTHoiwSRFTR8qEWq1ZKPmll/xun97InwutDzeMRe29ESCiSrpp6am4HA4\nUFpaCgDw+XxwuVzIyclBcXFx4B72RGR+USX9kSNHkJubG+hauN1uuFwudHR0oKioCG63W9cg5Vgs\nlqAHxZ8kSUEPrZZVKvSzkGifCz33zQzZmv4///kPTp06hSeeeAI//elPAQCNjY1obW0FAFRXV8Pp\ndIZN/Nra2sDPTqcTTqdTm6iJKMDr9cLr9Ua9vOy89w899BD279+Pa9eu4fDhw2hqakJGRgb8fj+A\nT7+ZMjMzA78HVhzHee8TcdAF3aDlgbxkn28/ms+6XO5F7N6fPHkSVqsVDodjzpUkUteJiGS69+fP\nn0djYyNOnTqF8fFxXLt2DVVVVbDZbBgYGEBWVhb6+/thtVrjFS+A+JzWEIUZpqKK19RbyfC50H0Y\n7sGDB9Hb24vOzk4cO3YMW7ZswQsvvICysjJ4PB4AgMfjQXl5uepAiCg+FJ2nn/mm3Lt3L1paWpCT\nk4MzZ85g7969ugRHRNpLyBtYsnuvHTN070PpObNutOtNZHK5l5DDcJnksYu0r8yyH7VMcjXDf5MV\nh+ESCYZJTyQYJj2RYBKypg8lSi2mhWS7bZeWtyWLddloltcKL60lIsWY9ESCSYruPcXOjF36SF1Y\nJXe4DV1e7hoRM87Co8ddldnSEwmGSU8kGCY9kWBY0wsmEYaeqqmttZrbQcv6Xw09tsOWnkgwTHoi\nwTDpiQTDmj4B6XW9uRnre6X0GoZr1styOQyXiGQx6YkEw+59AlDTRRd9lqFYSx+lM/QkUpnElp5I\nMEx6IsEw6YkEw5o+AWhZl5u93jSSklN2ak6TKtmuHF5aS0SymPREgmHSEwmGNX0CYl0evVjHMMTz\ndllKllczHdgMtvREgpFt6e12OxYuXIhbbrkFKSkpaGtrg8/nw9e+9jV0d3fDbrejoaEBixcvjke8\nRKSSbEtvsVjg9XrR3t6OtrY2AIDb7YbL5UJHRweKiorgdrt1D5SItBFVTR9aczQ2NqK1tRUAUF1d\nDafTGTbxa2trAz87nU44nc7YIxWYaOPltRTrGAatpt2KZl1qx1l4vV54vd7o45G7P/3nPvc5LFq0\nCLfccgu+9a1v4Rvf+AYyMjLg9/sDQWRmZgZ+D6xYx/vTi4ZJH39aHsiL13x70Q7UkW3pz507h+zs\nbHzwwQdwuVxYvXr1TRvS8luRiPQlm/TZ2dkAgNtuuw3bt29HW1sbbDYbBgYGkJWVhf7+flitVt0D\nTXR6zXYjty6KjZ7DcPW6LHfmtXKNcMQDeWNjYxgeHgYAjI6O4vTp08jPz0dZWRk8Hg8AwOPxoLy8\nPOZAiSi+Itb0nZ2d2L59OwBgcnISjzzyCPbt2wefz4fKykr09PTMecqONX0wtvSJTc//X7TrjZZc\n7skeyNNrw6Jh0ie2ZEp6DsPViZrENUuS88vlBi1n0tXrrjzR4jBcIsEw6YkEw6QnEgxrep0kwxRX\nWt0tJprlk4lRB2J5aS0RhcWkJxIMk55IMKzp4yQRa9pEPA6hF7l9EemyXD0H64Q7zqJq7D0RJR8m\nPZFg2L0nTeh5is4Md4Q1S2mjxZBdtvREgmHSEwmGSU8kGNb0NCc1daxWp6GUxmHEetXS6tp8DsMl\norCY9ESCYdITCYY1PelCi6mcZ2g9TFULRh134F1riUgxJj2RYJj0RIJhTU+mp9XxATU3kozXtQVK\nt8tLa4lIlu5Jr+S+2fHCmKLDmKJn1rjCYdKbRLLFNHMLc61vZa4mJkmSIj6UvDb0/W3evDnq9xr6\nWiUxRXpttNtn955IMEx6IsHoetdaIjKGIXetNctli0QUjN17IsEw6YkEw6QnEgyTnkgwTHoiwTDp\niQTzP9F1dFRExrmoAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106660d90>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAEHCAYAAABlS0A3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG0xJREFUeJzt3W9sE+cdB/CvodEgWwCn1E4CpdlaEpYQErcpaKhshswJ\nXUkIQ82otigCNE3T0EalCVHQSsY2Ela6Cra+mtjmdVpoyosSKE3JCI7aMha1C2KMStkoKSm1s9I4\nkD8wCLm96OLGxrnz+f46z/cjWYL4/NzPZ//8PM/d8zznkCRJAhEJY5rVARCRuZj0RIJh0hMJhklP\nJBgmPZFgmPREgmHSTwHTpk3D+++/b3UYqni9Xhw8eNDqMITEpLdIbm4u0tPTkZGRgaysLGzcuBHD\nw8OKr9OaLPX19aitrU369Xrt0+FwwOFwqC5r9+7dmDZtGtrb2/UKTzhMeos4HA4cO3YMg4OD+Pvf\n/4533nkHP//5zxN6ndb9pqqLFy/i8OHDyMnJsTqUlMakt4GcnBysXr0a58+fx8DAANasWQOXy4XM\nzExUVlbiypUrAICdO3fizTffxJYtW5CRkYEf/vCHkTLa2tqQl5cHp9OJLVu2TLovuQGYZ86cwfLl\ny+F0OlFSUoKOjo7Ic16vF88++ywee+wxzJo1CxUVFfjkk08iz//xj3/EAw88gLlz5+JnP/sZcnNz\ncfLkSbS2tqKhoQEvv/wyMjIy4PF4Iq/p6emZtLx4tmzZgr179yItLU12O5LHpLfQeAL29vbi9ddf\nx8MPP4yxsTFs3rwZly9fxuXLlzFz5sxIEv/iF7/AihUr8OKLL2JwcBAHDhyIlPXaa6/hnXfewblz\n59Dc3Iw33nhDVSxXrlzBmjVr8OyzzyIcDmPfvn1Yv359VCI2NTXhD3/4A/7zn//g1q1b2LdvHwDg\nwoUL+MEPfoCmpiYEg0Fcv34dH330ERwOB1avXo0dO3Zgw4YNGBwcRFdXV+S9//nPf45bXjyvvPIK\nZsyYgccff1zV+6K7MektIkkSqqur4XQ6sWLFCni9XuzYsQOZmZlYt24dZsyYgS984QvYsWNHVI07\n/tpY27dvx6xZs3D//fdj5cqVOHv2rKp4/vSnP+Eb3/gGVq9eDQD4+te/jtLSUrz22msAPu0WbNy4\nEQ899BBmzJiBmpqayD4OHz6MqqoqLF++HGlpadi9e3dUN0KSpLtidjgc2LRpU9zyYg0ODmLnzp3Y\nv3+/qvdE8d1jdQCicjgcOHLkCFatWhX195GRETz99NN44403EA6HAQBDQ0OQJCmSSPH65VlZWZF/\np6enY2hoSFU8H3zwAV555RUcPXo08rfR0dGo+CbuY+bMmZF9fPTRR5g/f37Uc/fee6/iPicrL9b4\nicAFCxZE/sZ5YsljTW8zzz//PLq7u9HZ2Ylr166ho6MjqqY06kTeggULUFtbi3A4HHkMDg5i27Zt\nimXm5OTgww8/jPz/xo0bUd0CrTG3t7fjwIEDyM7ORnZ2Nnp7e1FTU4PnnntOU7miYk1vM0NDQ5g5\ncyZmz56N/v5+/PSnP4163u124+LFi7JlyNWCkiRhbGwM//3vf6N+SL7zne/g0UcfxYkTJ1BWVobb\nt2/jzJkzWLhwIebNmydb7vr16/GVr3wFf/3rX/HII4+gvr4+atusrCz85S9/iWqtKMU50cmTJzE6\nOhp5zaOPPooXXngh0hUhdVjT28zWrVtx48YNzJ07F8uXL8fjjz8elSg/+tGPcPjwYWRmZmLr1q1x\ny5C7Bu5wONDU1ISZM2ciPT0d6enpWLhwIebPn48jR45gz549cLlcWLBgAZ5//vmoxJxY5sR9FBYW\n4te//jU2bNiAnJwcZGRkwOVy4XOf+xwA4MknnwQA3HvvvSgtLVUsL1ZmZiZcLhdcLhfcbjemT58O\np9OJz3/+87LHkuJzcBEN0tvQ0BCcTif+/e9/44EHHrA6HIrBmp50cfToUYyMjGB4eBg//vGPsWTJ\nEia8TTHpSRctLS2YN28e5s2bh4sXL+LQoUNWh0STYPOeSDCG1vStra1YtGgRFi5ciL179xq5q0lt\n2rQJbrcbRUVFkb/19/fD5/MhLy8P5eXlGBgYMDWm3t5erFy5EoWFhVi8eHFkZJ2Vcd28eRPLli1D\nSUkJCgoK8Mwzz1ge07g7d+7A4/GgsrLSNjHl5uZiyZIl8Hg8WLp0qW3iSoRhSX/nzh1s2bIFra2t\nuHDhApqamvDee+8ZtbtJbdy4Ea2trVF/a2xshM/nQ3d3N8rKytDY2GhqTGlpaXjhhRfwz3/+E2fO\nnMGLL76I9957z9K4ZsyYgVOnTuHs2bM4d+4cTp06hbfeesvyYwUA+/fvR0FBQeTsvh1icjgcCAQC\n6OrqQmdnp23iSohkkNOnT0sVFRWR/zc0NEgNDQ1G7U7WpUuXpMWLF0f+n5+fL4VCIUmSJCkYDEr5\n+fmWxDVu7dq1Ultbm23iGh4elkpLS6Xz589bHlNvb69UVlYmtbe3S2vWrJEkyR6fX25urnT16tWo\nv9khrkQYVtNfuXIF999/f+T/8+fPj8wWs1pfXx/cbjeATwe79PX1WRZLT08Purq6sGzZMsvjGhsb\nQ0lJCdxud6T7YXVMTz/9NJ577jlMm/bZV9XqmIBPa/rx+Qm//e1vbRNXIgwbkZcq87aTXcxBD0ND\nQ1i/fj3279+PjIwMy+OaNm0azp49i2vXrqGiogKnTp2yNKZjx47B5XLB4/EgEAjE3caqz+/tt99G\ndnY2Pv74Y/h8PixatMgWcSXCsJp+3rx56O3tjfy/t7c3alKGldxuN0KhEAAgGAzC5XKZHsPt27ex\nfv161NbWorq62jZxAcDs2bPxxBNP4N1337U0ptOnT6OlpQVf/OIX8dRTT6G9vR21tbW2OE7Z2dkA\ngPvuuw/r1q1DZ2enLeJKhGFJX1pain/961/o6enBrVu38PLLL6Oqqsqo3alSVVUFv98PAPD7/ZGk\nM4skSdi8eTMKCgqihtJaGdfVq1cjZ5tv3LiBtrY2eDweS2Pas2cPent7cenSJRw6dAirVq3CSy+9\nZPnnNzIygsHBQQDA8PAwTpw4gaKiIsvjSpiRJwyOHz8u5eXlSQ8++KC0Z88eI3c1qQ0bNkjZ2dlS\nWlqaNH/+fOl3v/ud9Mknn0hlZWXSwoULJZ/PJ4XDYVNjevPNNyWHwyEVFxdLJSUlUklJifT6669b\nGte5c+ckj8cjFRcXS0VFRdIvf/lLSZIky4/VuEAgIFVWVtoipvfff18qLi6WiouLpcLCwsh32+q4\nEsXBOUSC4TBcIsEknfR2GG1HRElIpk8wOjoqPfjgg9KlS5ekW7duScXFxdKFCxeitgHABx98WPSQ\nk1RN39nZiYceegi5ublIS0vDhg0bcOTIkWSKSpj0/yWj9D4FMbFcpbK1bKvmtXqWpWXbXbt2Ke5b\nj4ea9xsbk9pjZdTDrGMV73Hq1Cns2rUr8lCS1OCceKPt/va3vyVTFBFp5PV64fV6I/+PXWItVlI1\nvV1HGhGRsqSS3orRduPDGh0Oh+qmcrzm0Pi/J5ar9GOmtO3EfcRuG/vQ0txXQ2m/cjFPrD3UHgst\nMcmJjcms46hE7ljZTVLX6UdHR5Gfn4+TJ08iJycHS5cuRVNTE7785S9/VrCBrYHYkJX2JfcW9Yxz\n4n7UxiS3vdJHpOY9qDkWem2r52tJ2fgP6WSS6tPfc889+M1vfoOKigrcuXMHmzdvjkp4IrIvw0bk\nmdnv11ITWlXjytHSktGz1SMXB2vnxJndslGq6Tkij0gwTHoiwTDpiQQj/L3stJyxjmVU31qJWecs\nROrHT+WrE6zpiQTDpCcSDJOeSDAp2ac38jq2FmquY+vVD1fal5r3a/e+aCL06otree92P26s6YkE\nw6QnEgyTnkgwhvbp1fZNE6XlOrbamPSawaa0rZbr51qYdb5jIrX7MStGu/fF9cKankgwTHoiwRja\nvJe7fKLnIhNaYlJ63qwmn5b9aDk2qTZdVul7oaXblGrHIlms6YkEw6QnEgyTnkgwpg7D1WuYKkUz\najqwnv1aNZdN9eqXq4kpXll2pEfMrOmJBMOkJxIMk55IMCk5tVYLI/uTdqC2z2fU+Ae9ltqONdU/\nv1hGnHdgTU8kGCY9kWCY9ESCMbVPb8ex93ruR5Sx24C2vqaWfrld504YxYj4WdMTCUY26Tdt2gS3\n242ioqLI3/r7++Hz+ZCXl4fy8nIMDAwYHiQR6Uc26Tdu3IjW1taovzU2NsLn86G7uxtlZWVobGw0\nJDCHwxH1kCQp8lAycdtELllNfBgldj9yDzPFHiu54yYXo5r3oPb9ysWo9KC7ySb9ihUr4HQ6o/7W\n0tKCuro6AEBdXR1effVV46IjIt2pPpHX19cHt9sNAHC73ejr69M9KCJKXCAQQCAQSHh7h6TQBurp\n6UFlZSX+8Y9/AACcTifC4XDk+czMTPT3999dcJxmm5YFDu1w5l/NfrTEoOcqOkbNlFNi1fsnRLrD\nk1Fd07vdboRCIWRlZSEYDMLlck26rZ5JodeUST1/AOTKmgrLbukZs5opvHoOw6W7qb5kV1VVBb/f\nDwDw+/2orq7WPSgiMo5s8/6pp55CR0cHrl69Crfbjd27d2Pt2rWoqanB5cuXkZubi+bmZsyZM+fu\nguM0MfRq/qplxdrucjHYlZ4Dl+wy4UZESs17xT69njtm0tsbk35q0L1Pr3bniT6X7N1F4zHzB2Uq\nUZNQXOLKHJxaS0SaMemJBMOkJxJMyiyXZcfls/U8UZkKfVO93t9UX+JKiZZzJXpgTU8kGCY9kWAM\nbd7b5Tq9UYwa7mvXy1Jamuh2HLJsFjXvz4z3zpqeSDBMeiLBMOmJBGPZMFwjGTWkV8/pv1pYNTed\nc+ITY/cVe1nTEwmGSU8kGCY9kWBsc4cbI+dJG3XXGtGW3jJrGC4ZizU9kWCY9ESCsWwYbiyzVmJV\n2q/ca1NxqLAWeg4zVjNLMhWp6fpYPauQNT2RYJj0RIJh0hMJJmVWztGL2nMHVqzCY5c+oGgrA8lR\n872x6vNL9LvKmp5IMEx6IsEw6YkEY+rUWquG4aYCo4b7GnUnoan++Wh5f3pep1djvFyl8ljTEwlG\nNul7e3uxcuVKFBYWYvHixThw4AAAoL+/Hz6fD3l5eSgvL8fAwIApwRKRdrJ3rQ2FQgiFQigpKcHQ\n0BAeeeQRvPrqq/j973+PuXPnYtu2bdi7dy/C4TAaGxujC47TxLBD896ul+wm0vNusUbeKDTZclOB\nlu+Jlua9HsdR6a61kFRYu3at1NbWJuXn50uhUEiSJEkKBoNSfn7+XdsCUPVQ83ot9CxLqexk35+W\nh9r3a9RxTkVGfV/NPq5K+0j4RF5PTw+6urqwbNky9PX1we12AwDcbjf6+voSLYaIdBYIBBAIBBLe\nXrZ5P25oaAhf+9rX8JOf/ATV1dVwOp0Ih8OR5zMzM9Hf3x9dsMpmcWwYZo1a0lKWUtly+zGq26C0\nH7OOcyrSsxul5rV6U2reK9b0t2/fxvr161FbW4vq6moAn9buoVAIWVlZCAaDcLlcqgNT++XUi54/\nGFpea9SPgJErsU716bFaTDwedv+hlT17L0kSNm/ejIKCAmzdujXy96qqKvj9fgCA3++P/BgQkf3J\nNu/feustfPWrX8WSJUsiv04NDQ1YunQpampqcPnyZeTm5qK5uRlz5syJLljlBIxUGOyhZ5PcrOa+\n2jjkTPWaXq8mutU1vVLzPqE+fbI7lsOkZ9LbjShJb+pyWRPZ/S4g8SjFYIcTeXqy42egJz2HiSe7\nrZ44tZaI4mLSEwmGSU8kGOGWy9KTliWUjMLBOdHUnHw06w4+Rt0ZaXw/nFpLRFGY9ESCYdITCcbU\n5bIm0jI4Z6pfP9bCzP6kHem1zJXa42jHW7ZNhjU9kWCY9ESCseySnZHNTDvNbU5mP1rnLagpK9Wp\nGc6tZRh1rFTuUrKmJxIMk55IMEx6IsGYOrXWqOmIVl3Cs+PSW1O9D68nue+ckcNwrcaankgwTHoi\nwTDpiQRj2V1rlbad6rQcC9GG1k6k5XZoWvrlauIysw+fzH5Z0xMJhklPJBgmPZFgUma5rFS6bVAy\njBojnor0nGad7NgQLdfpjaTHGADW9ESCYdITCSZlhuFqaXrJsfuQSRHotdqNUll2We1GCz3iYE1P\nJBjZpL958yaWLVuGkpISFBQU4JlnngEA9Pf3w+fzIS8vD+Xl5RgYGDAlWCLSTvGutSMjI0hPT8fo\n6Cgee+wx7Nu3Dy0tLZg7dy62bduGvXv3IhwOo7GxMbrgOHfO1GuWnZbmvV3vlmvHs/d2PCNt5Ig8\nOXZp3idC6a61kBI0PDwslZaWSufPn5fy8/OlUCgkSZIkBYNBKT8//67tVRRtKwCSfsiVo7QfNTGZ\nRcuxSPa4KcVhFbWfp5r3q/e2SsdJ8UTe2NgYHn74YVy8eBHf//73UVhYiL6+PrjdbgCA2+1GX19f\n3NfW19dH/u31euH1epV2R0RJmphvchSb9+OuXbuGiooKNDQ04Jvf/CbC4XDkuczMTPT390cXrNTE\nsCm9moBamqFKMZl1XO3SjbBqMstkMcSLQ88uptZtlXIv4bP3s2fPxhNPPIF3330XbrcboVAIABAM\nBuFyuRIthogsJpv0V69ejZyZv3HjBtra2uDxeFBVVQW/3w8A8Pv9qK6uVr1jh8MR9TBL7H5jH5Ik\nRT0mknsutuzYbdW8XzUxkTm0fAZK34VEv1PxvgvxylUi26cPBoOoq6vD2NgYxsbGUFtbi7KyMng8\nHtTU1ODgwYPIzc1Fc3OzqoNARNZJuE+vumCFfoVd+6lyfSij+mKx29vxMqKRUqFPr0Sv+zZqmdiT\n6HEydZadHVYXUbttsh/IVJwJqBctP4h2paaCU/u83jgMl0gwTHoiwTDpiQRjap/eqpN1Wk6U6DUm\n3i4n56yi5jOIxfMfn0l0cI4c1vREgmHSEwmGSU8kmJRZDVcNuX4e++XWsMNKs3p+fmat4aBljMlk\nWNMTCYZJTyQYJj2RYKZEn97Itc+MGns/1ajpp+t5nV4pDqNY9XlqWbxjHGt6IsEw6YkEkzJTa42a\nr6zn2mepMO9bL2ovQ5k9fdTO9JyGO9lqO3JY0xMJhklPJBgmPZFgDO3Tm3XZSs+lqYxavywV6Dk8\n1Kr19USehpso1vREgmHSEwmGSU8kGEP79FYNiTSqbzoV+4Sp9v70GkeRSNlWMGP8A2t6IsEw6YkE\nw6QnEkzKTK1lv9wYyR5Xq8be2+V2YWYtvWXEcWRNTySYhJL+zp078Hg8qKysBAD09/fD5/MhLy8P\n5eXlkXvYE5H9JZT0+/fvR0FBQaRp0djYCJ/Ph+7ubpSVlaGxsdHQIIFPmzXjDyWSJEU9tGybaDmp\nYuJxVHrEHhuzjoWafcXGrGZbuYeZtBzXZF6r2Kf/8MMPcfz4cezcuRO/+tWvAAAtLS3o6OgAANTV\n1cHr9cZN/Pr6+si/vV4vvF5vwoERUWICgQACgUDC2zskhZ+IJ598Ejt27MD169exb98+HD16FE6n\nE+FwGMCnvzSZmZmR/0cK/n9NoRc166breT+6qVLDjzNqvTk9T+TZYQEOq74nao/jZBPC5OKRbd4f\nO3YMLpcLHo9n0kKsaA4RUfJkm/enT59GS0sLjh8/jps3b+L69euora2F2+1GKBRCVlYWgsEgXC5X\n3Nfb5ddQS1mpxsi78qopV8trtUxvTnVqv8vxPj9Ny2Xt2bMHvb29uHTpEg4dOoRVq1bhpZdeQlVV\nFfx+PwDA7/ejurpadidEZB+qrtOP/4Js374dbW1tyMvLQ3t7O7Zv325IcESkP8UTeUkXHKeJYYdR\nS1pm5KUCq5r3SsxavUgvdjnhm8wJUKWT6LadWmvWENBUT3JAv3MYRn6RU23prVSYhpssDsMlEgyT\nnkgwTHoiwdh2aq3Rt/5JZXou0y3ytGOjlt7S81yWESc5WdMTCYZJTyQY2zbvaXJyzUe115flnjPq\nkmssq67T2+VuRlo+E66GS0SKmPREgmHSEwkmJfv0U338fCwt/dhElgvTYz9q92sWNe9Pz7vlTBaD\nkXSZWktEUw+TnkgwTHoiwdi2T2/VFFCrGDXf3Kw7z+jJqKW3psJ7V3NNfzKs6YkEw6QnEoxtmvdq\nmuxTvTmvtK1Ry4PZsbkby4zmr1ZquxFmf59Z0xMJhklPJBgmPZFgbNOnnwr9dDlm9cuV9iv3vGjD\nm/Wk5dyIXtOQOQyXiOJi0hMJhklPJBjb9OmnOjVTXGNp6ZdrYeSqrlrKkmPVMFwt+9Ey/Zer4RKR\nIsWaPjc3F7NmzcL06dORlpaGzs5O9Pf341vf+hY++OAD5Obmorm5GXPmzDEjXiLSSLGmdzgcCAQC\n6OrqQmdnJwCgsbERPp8P3d3dKCsrQ2Njo+GBEpE+EurTx/ZRWlpa0NHRAQCoq6uD1+uNm/j19fWR\nf3u9Xni93uQjTTFmLXGltC+jrvGbec3eDuMD1BxzLdfpkznvEAgEEAgEEt+f0v3pv/SlL2H27NmY\nPn06vve97+G73/0unE4nwuFwJMjMzMzI/ycGa4cPyyp6Jr2agT163lddS9LbZWCPXifyzJqLr9dJ\nQE33p3/77beRnZ2Njz/+GD6fD4sWLborqFSYnUVEn1JM+uzsbADAfffdh3Xr1qGzsxNutxuhUAhZ\nWVkIBoNwuVyaA7FLzaCGns3qic/rudpNKh5XPSX7frV8flqo/bx0v8PNyMgIBgcHAQDDw8M4ceIE\nioqKUFVVBb/fDwDw+/2orq5WvWMisoZsTd/X14d169YBAEZHR/Htb38b5eXlKC0tRU1NDQ4ePBi5\nZEdEqUHxRF7SBas8kZeKzVCjmvdqy9KyXzlT4USeHLM+Py2Sbd5rOpFnlFT4UsSy4xJXasvSwqq7\nuhpFr6HPqYbDcIkEw6QnEgyTnkgwlvXpU6HPByQ/osusJa6UytJCz/EBcuXSZ4ycljuONT2RYJj0\nRIJh0hMJhstlKUj29lp26ZdroWVpJj3fj17Lgek5q86ssfdqtucS2EQUF5OeSDBs3ivQq2mp5rV2\nbOoD5jVh1XSb9Jx6qtdcEbXvz+yuAWt6IsEw6YkEw6QnEsyU7NPrOYdcr+WWUpFRc/6V6DlkWcvd\nZNXQsjip2VjTEwmGSU8kGCY9kWCmRJ9ezyGtU6EvroaaqcNm9a2NOj+gpa9t1pp4ZmBNTyQYJj2R\nYJj0RIKZEn16kfvhgDXjENSyqm+t13V6o85vaI0rGazpiQRjeNKruW+2WRhT6rLrcbJrXPEw6U0y\nfkvvyW7tLRdT7GslSYp6qKHmtbt27YraNjYOvWJSY+XKlbL7jX3IHXM1McuV63A4ouJSInfszMDm\nPZFgmPREgjH0rrVEZA1L7lor2mU0olTB5j2RYJj0RIJh0hMJhklPJBgmPZFgmPREgvkfw2Lxr+kk\nCg8AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1080c35d0>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAEHCAYAAABlS0A3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFFVJREFUeJzt3XtQVOUfx/HPUZmUQFzAPYBUWyoaiLDl6IxdPEYLlUGY\nSTZFpE41zThlM2ZoU2yXkbVf2djlj6bpsuWMZU2jYEaSskxpxVSaXYcyCVSWSXdBriHw/P5wPLWB\nu8tlOWvP5zXDDLvsPufr5c05Z108ihBCgIikMcboAYhodDF6IskweiLJMHoiyTB6IskweiLJMPrz\n0JgxY/D7778bPcagaJqG119/3egxCIx+1FgsFkRGRiI6OhoJCQlYvnw52tvbAz5vuLHY7XYUFhYO\n+fkjtU1FUaAoSlDPP336NG677TZceumlGDNmDKqrq/s95tFHH0V8fDzi4+NRXFw8InPLgtGPEkVR\nsHPnTrS2tuLbb7/F119/jWeeeSao5w13u+eja6+9Flu2bEFCQkK/X8Orr76KHTt24NChQzh06BDK\ny8vx6quvGjTp+YfRGyApKQk33HADfvjhBzQ3N+Pmm2+G2WxGbGwscnNzcezYMQDAY489hs8++wyr\nVq1CdHQ0HnzwQX2NyspKpKSkwGQyYdWqVefclr83XH755ZeYP38+TCYTMjMzffaomqbhiSeewNVX\nX42JEyciJycHJ0+e1L/+9ttv45JLLkF8fDyefvppWCwW7NmzBxUVFSgtLcV7772H6OhoWK1W/Tl1\ndXXnXO+fIiIi8OCDD+Kqq67C2LFj+33d6XRizZo1SEpKQlJSEtasWYO33nrrnL9O+hdBo8JisYhP\nP/1UCCFEfX29SEtLE0888YQ4efKk+PDDD0VnZ6dobW0VS5cuFfn5+frzNE0Tr7/+us9aiqKI3Nxc\n0dLSIurr68XkyZNFRUXFgNstKSkRd911V7/7jx49KuLi4sTHH38shBCisrJSxMXFiRMnTgghhFiw\nYIGYNm2a+PXXX0VnZ6fQNE0UFxcLIYT48ccfRVRUlNi3b5/o7u4Wa9asEREREWLPnj1CCCHsdrso\nLCz02d6CBQvE1KlTB1zPn+TkZFFdXe1zX0xMjKipqdFvf/311yI6OjrgWnQG9/SjRAiB/Px8mEwm\nXHPNNdA0DevXr0dsbCwWL16M8ePHIyoqCuvXr+93DisG2FsXFxdj4sSJuOiii7Bw4UIcPHhwUPNs\n2bIFN910E2644QYAwPXXX485c+bgo48+AnDmtGD58uWYNm0axo8fj4KCAn0bH3zwAfLy8jB//nxE\nRETgqaee8jkEF0L0m1lRFKxYsWLA9Qarra0NMTEx+u2JEyeira1tSGvJaJzRA8hCURTs2LED1113\nnc/9HR0dePjhh/HJJ5/A6/UCOPOXWgihhzTQeXlCQoL+eWRk5KD/0v/xxx94//33UV5ert/X09Pj\nM98/tzFhwgR9G8ePH0dycrLP1+Li4gJu81zrDVZUVBROnTql325paUFUVNSQ1pIR9/QGe/7551Fb\nW4uamhq0tLSgurraZ08ZqhfyLr74YhQWFsLr9eofra2tWLt2bcA1k5KScPToUf12Z2enz/l5qF88\nTEtL8zlK+O677zBr1qyQbvO/hNEbrK2tDRMmTEBMTAw8Hg+efPJJn6+rqorDhw/7XWOgw/9/fq2v\nrw9//fUXurq60NXVhb/++gt33XUXysvLsXv3bvT29qKrqwsul0t/EdHfukuWLEF5eTm++OILdHd3\nw263+zw2ISEBdXV1/Z7vb85/Ozvvvz8HgLvvvhubNm3C8ePHcezYMWzatAn33HNP0GvLjtEbbPXq\n1ejs7ER8fDzmz5+PG2+80WdP+dBDD+GDDz5AbGwsVq9ePeAa/v4NXFEUbN26FRMmTEBkZCQiIyMx\nffp0JCcnY8eOHdiwYQPMZjMuvvhiPP/88z5h/nPNf24jLS0NL730EpYtW4akpCRER0fDbDbjggsu\nAAAsXboUABAXF4c5c+YEXG8gM2bMQGRkJI4fP46cnBxceOGFqK+vBwDcf//9yM3NRXp6OmbPno3c\n3Fzcd99951yLfCliMN9+iQbQ1tYGk8mE3377DZdcconR41AA3NPTkJSXl6OjowPt7e1Ys2YNZs+e\nzeDPE4yehqSsrAxTpkzBlClTcPjwYbz77rtGj0RB4uE9kWRCuqevqKjAzJkzMX36dGzcuDGUmzqn\nFStWQFVVpKen6/d5PB7YbDakpKQgOzsbzc3NozpTQ0MDFi5ciLS0NMyaNQsvvvii4XN1dXVh3rx5\nyMzMRGpqKtatW2f4TGf19vbCarUiNzc3bGayWCyYPXs2rFYr5s6dGzZzBSNk0ff29mLVqlWoqKjA\nTz/9hK1bt+Lnn38O1ebOafny5aioqPC5z+FwwGazoba2FllZWXA4HKM6U0REBF544QX8+OOP+PLL\nL/HKK6/g559/NnSu8ePHo6qqCgcPHsShQ4dQVVWFzz//3PDfKwDYvHkzUlNT9Vf7w2EmRVHgcrlw\n4MAB1NTUhM1cQQnV+3v3798vcnJy9NulpaWitLQ0VJvz68iRI2LWrFn67RkzZgi32y2EEKKxsVHM\nmDHDkLnOuuWWW0RlZWXYzNXe3i7mzJkjfvjhB8NnamhoEFlZWWLv3r3i5ptvFkKEx5+fxWLRf07h\nrHCYKxgh29MfO3YMF110kX47OTnZ540fRmpqaoKqqgDOvPmlqanJsFnq6upw4MABzJs3z/C5+vr6\nkJmZCVVV9dMPo2d6+OGH8b///Q9jxvz9V9XomYAze/qzP6/w2muvhc1cwQjZe+/Pl5/jHsx/7jDS\n2trasGTJEmzevBnR0dGGzzVmzBgcPHgQLS0tyMnJQVVVlaEz7dy5E2azGVarFS6Xa8DHGPXnt2/f\nPiQmJuLPP/+EzWbDzJkzw2KuYIRsTz9lyhQ0NDTotxsaGnx+SMNIqqrC7XYDABobG2E2m0d9htOn\nT2PJkiUoLCxEfn5+2MwFADExMVi0aBG++eYbQ2fav38/ysrKcOmll+KOO+7A3r17UVhYGBa/T4mJ\niQCAyZMnY/HixaipqQmLuYIRsujnzJmDX3/9FXV1deju7sZ7772HvLy8UG1uUPLy8uB0OgGc+Q8Z\nzkY3WoQQWLlyJVJTU33eWmvkXCdOnNBfbe7s7ERlZSWsVquhM23YsAENDQ04cuQI3n33XVx33XV4\n5513DP/z6+joQGtrKwCgvb0du3fvRnp6uuFzBS2ULxjs2rVLpKSkiKlTp4oNGzaEclPntGzZMpGY\nmCgiIiJEcnKyeOONN8TJkydFVlaWmD59urDZbMLr9Y7qTJ999plQFEVkZGSIzMxMkZmZKT7++GND\n5zp06JCwWq0iIyNDpKeni2effVYIIQz/vTrL5XKJ3NzcsJjp999/FxkZGSIjI0OkpaXpf7eNnitY\nfHMOkWT4NlwiyQw5+nB4tx0RDcFQzgl6enrE1KlTxZEjR0R3d7fIyMgQP/30k89jAPCDH/ww6MOf\nIe3pa2pqMG3aNFgsFkRERGDZsmXYsWNHv8eVlJRgwYIFKCkpQVVVlf7fQBn9UVJSYvgMnOm/M5PR\nc1VVVaGkpET/CGRIb84Z6N12X331Vb/H2e12/YOIQkPTNGiapt/+93+59m9D2tOH6zuNiCiwIUU/\nmHfb/fM7ULjgTMHhTMEL17kGMqR/p+/p6cGMGTOwZ88eJCUlYe7cudi6dSsuv/zyvxdWFAxhaSIa\npkDtDemcfty4cXj55ZeRk5OD3t5erFy50id4IgpfIXtHHvf0RMYI1B7fkUckGUZPJBlGTyQZRk8k\nGUZPJBlGTyQZRk8kGUZPJBlGTyQZRk8kGUZPJBlGTyQZRk8kGUZPJBlGTyQZRk8kGUZPJBlGTyQZ\nRk8kGUZPJBlGTyQZRk8kGUZPJBlGTyQZRk8kGUZPJBlGTyQZRk8kGb/Rr1ixAqqqIj09Xb/P4/HA\nZrMhJSUF2dnZaG5uDvmQRDRy/Ea/fPlyVFRU+NzncDhgs9lQW1uLrKwsOByOkA5IRCMr4KWq6+rq\nkJubi++//x4AMHPmTFRXV0NVVbjdbmiahl9++aX/wrxUNZEhArU3brALNjU1QVVVAICqqmhqajrn\nY+12u/65pmnQNG2wmyOiAFwuF1wuV9CPH/Se3mQywev16l+PjY2Fx+PpvzD39ESGCNTeoF+9P3tY\nDwCNjY0wm81Dn46IRt2go8/Ly4PT6QQAOJ1O5Ofnj/hQRBQ6fg/v77jjDlRXV+PEiRNQVRVPPfUU\nbrnlFhQUFKC+vh4WiwXbtm3DpEmT+i/Mw3siQwRqL+A5fag2TEShMeLn9ER0fmP0RJJh9ESSYfRE\nkmH0RJJh9ESSYfREkmH0RJJh9ESSYfREkmH0RJJh9ESSYfREkmH0RJJh9ESSYfREkmH0RJJh9ESS\nYfREkmH0RJJh9ESSYfREkmH0RJJh9ESSYfREkmH0RJJh9ESSYfREkvEbfUNDAxYuXIi0tDTMmjUL\nL774IgDA4/HAZrMhJSUF2dnZaG5uHpVhiWj4/F611u12w+12IzMzE21tbbjyyiuxfft2vPnmm4iP\nj8fatWuxceNGeL1eOBwO34V51VoiQwzrqrUJCQnIzMwEAERFReHyyy/HsWPHUFZWhqKiIgBAUVER\ntm/fPoIjE1EojQv2gXV1dThw4ADmzZuHpqYmqKoKAFBVFU1NTQM+x263659rmgZN04Y1LBH153K5\n4HK5gn6838P7s9ra2rBgwQI8/vjjyM/Ph8lkgtfr1b8eGxsLj8fjuzAP74kMMazDewA4ffo0lixZ\ngsLCQuTn5wM4s3d3u90AgMbGRpjN5hEal4hCzW/0QgisXLkSqampWL16tX5/Xl4enE4nAMDpdOrf\nDIgo/Pk9vP/8889x7bXXYvbs2VAUBQBQWlqKuXPnoqCgAPX19bBYLNi2bRsmTZrkuzAP74kMEai9\noM7pQ7FhIgqNYZ/TE9F/C6MnkgyjJ5IMoyeSDKMnkgyjJ5IMoyeSDKMnkgyjJ5IMoyeSDKMnkgyj\nJ5IMoyeSDKMnkgyjJ5IMoyeSDKMnkgyjJ5IMoyeSDKMnkgyjJ5IMoyeSDKMnkgyjJ5IMoyeSDKMn\nkgyjJ5IMoyeSjN/ou7q6MG/ePGRmZiI1NRXr1q0DAHg8HthsNqSkpCA7OxvNzc2jMiwRDV/Aq9Z2\ndHQgMjISPT09uPrqq/Hcc8+hrKwM8fHxWLt2LTZu3Aiv1wuHw+G7MK9aS2SIYV+1NjIyEgDQ3d2N\n3t5emEwmlJWVoaioCABQVFSE7du3j9C4RBRq4wI9oK+vD1dccQUOHz6MBx54AGlpaWhqaoKqqgAA\nVVXR1NQ04HPtdrv+uaZp0DRtRIYmor+5XC64XK6gHx/w8P6slpYW5OTkoLS0FLfeeiu8Xq/+tdjY\nWHg8Ht+FeXhPZIhhH96fFRMTg0WLFuGbb76Bqqpwu90AgMbGRpjN5uFPSkSjwm/0J06c0F+Z7+zs\nRGVlJaxWK/Ly8uB0OgEATqcT+fn5oZ+UiEaE38P777//HkVFRejr60NfXx8KCwvxyCOPwOPxoKCg\nAPX19bBYLNi2bRsmTZrkuzAP74kMEai9oM/pR3rDRBQaI3ZOT0T/DYyeSDKMnkgyjJ5IMoyeSDKM\nnkgyjJ5IMoyeSDKMnkgyjJ5IMoyeSDKMnkgyjJ5IMoyeSDKMnkgyjJ5IMoyeSDKMnkgyjJ5IMoye\nSDKMnkgyjJ5IMoyeSDKMnkgyjJ5IMoyeSDKMnkgyjJ5IMoyeSDJBRd/b2wur1Yrc3FwAgMfjgc1m\nQ0pKCrKzs/Vr2BNR+Asq+s2bNyM1NRWKogAAHA4HbDYbamtrkZWVBYfDEdIhiWjkjAv0gKNHj2LX\nrl147LHHsGnTJgBAWVkZqqurAQBFRUXQNG3A8O12u/65pmnQNG1kpiYincvlgsvlCvrxivB39XoA\nS5cuxfr163Hq1Ck899xzKC8vh8lkgtfrBQAIIRAbG6vf1hdWFARYmohCIFB7fg/vd+7cCbPZDKvV\nes5FFEXRD/uJKPz5Pbzfv38/ysrKsGvXLnR1deHUqVMoLCyEqqpwu91ISEhAY2MjzGbzaM1LRMMU\n8PD+rOrqav3wfu3atYiLi8Ojjz4Kh8OB5ubmfuf0PLwnMsawDu8HWgwAiouLUVlZiZSUFOzduxfF\nxcXDm5KIRk3Qe/pBL8w9PZEhRnRPT0TnP0ZPJBlGTyQZRk8kGUZPJBlGTyQZRk8kGUZPJBlGTyQZ\nRk8kGUZPJBlGTyQZRk8kGUZPJBlGTyQZRk8kGUZPJBlGTyQZRk8kGUZPJBlGTyQZRk8kGUZPJBlG\nTyQZRk8kGUZPJBlGTyQZRk8kGb/XpwcAi8WCiRMnYuzYsYiIiEBNTQ08Hg9uv/12/PHHH7BYLNi2\nbRsmTZo0GvMS0TAF3NMrigKXy4UDBw6gpqYGAOBwOGCz2VBbW4usrKx+16YnovAVcE8PoN9lb8vK\nylBdXQ0AKCoqgqZpA4Zvt9v1zzVNg6ZpQ5+UiAbkcrngcrmCfnzA69NfdtlliImJwdixY3H//ffj\n3nvvhclkgtfrBXDmG0JsbKx+W1+Y16cnMkSg9gLu6fft24fExET8+eefsNlsmDlzZr8NKIoy/EmJ\naFQEPKdPTEwEAEyePBmLFy9GTU0NVFWF2+0GADQ2NsJsNod2SiIaMX6j7+joQGtrKwCgvb0du3fv\nRnp6OvLy8uB0OgEATqcT+fn5oZ+UiEaE33P6I0eOYPHixQCAnp4e3HnnnVi3bh08Hg8KCgpQX19/\nzn+y4zk9kTECtRfwhbxQbZiIQiNQe3xHHpFkGD2RZBg9kWQYPZFkGD2RZBg9kWQYPZFkGD2RZBg9\nkWQYPZFkGD2RZBg9kWQYPZFkGD2RZBg9kWQYPZFkGD2RZBg9kWQYPZFkGD2RZBg9kWQYPZFkGD2R\nZBg9kWQYPZFkQh79YK6bPVo4U3A4U/DCda6BMPowwZmCE44zAeE710B4eE8kGUZPJJmQXrWWiIzh\nL+txRmyUiIzDw3siyTB6IskweiLJMHoiyTB6IskweiLJ/B/eXzVB4LOfOwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1080c3590>"
       ]
      }
     ],
     "prompt_number": 40
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Lab 7"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "import numpy as np\n",
      "\n",
      "\n",
      "def rowswap(n, j, k):\n",
      "    \"\"\" Swaps two\n",
      "    INPUTS: n -> matrix size\n",
      "            j, k -> the two rows to swap\n",
      "    \"\"\"\n",
      "    out = np.eye(n)\n",
      "    out[j, j] = 0\n",
      "    out[k, k] = 0\n",
      "    out[j, k] = 1\n",
      "    out[k, j] = 1\n",
      "    return out\n",
      "\n",
      "\n",
      "def cmult(n, j, const):\n",
      "    \"\"\" Multiplies a row by a constant\n",
      "    INPUTS: n -> array size\n",
      "            j -> row\n",
      "        const -> constant\n",
      "    \"\"\"\n",
      "    out = eye(n)\n",
      "    out[j, j] = const\n",
      "    return out\n",
      "\n",
      "\n",
      "def cmultadd(n, j, k, const):\n",
      "    \"\"\"Multiplies a row (k) by a constant and adds the result to\n",
      "    \ufffc\ufffc\ufffc\ufffcanother row (j)\"\"\"\n",
      "    out = eye(n)\n",
      "    out[j, k] = const\n",
      "    return out\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 41
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def ref(matrix):\n",
      "    \"\"\"\n",
      "    Put an n x (n+1) matrix in REF.\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    matrix : array_like, dtype=float, shape=(n, n+1)\n",
      "        The matrix representing an augmented linear system.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    ref : array_like, dtype=float, shape=(n, n+1)\n",
      "        The REF of the input matrix\n",
      "\n",
      "    \"\"\"\n",
      "    n = matrix.shape[0]\n",
      "    for i in range(0, n):\n",
      "        for j in range(i, n):\n",
      "            matrix = np.dot(cmultadd(n, j, i, -matrix[j,i] / matrix[i,i]), matrix)\n",
      "    for i in range(0, n):\n",
      "        matrix[i,:] /= matrix[i,i]\n",
      "    return matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 42
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "A = np.array([[4, 5, 6, 3], [2, 4, 6, 4], [7, 8, 0, 5]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 43
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ref(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 44,
       "text": [
        "array([[ 1.        ,  1.25      ,  1.5       ,  0.75      ],\n",
        "       [-0.        ,  1.        ,  2.        ,  1.66666667],\n",
        "       [ 0.        ,  0.        ,  1.        , -0.11111111]])"
       ]
      }
     ],
     "prompt_number": 44
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 2"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def LUDecomp(mat):\n",
      "    \"\"\"\n",
      "    Perform LU decomposition of a matrix using only elementary matrices\n",
      "    \"\"\"\n",
      "    n = mat.shape[0]\n",
      "    EL = []\n",
      "    L = np.eye(n)\n",
      "    U = mat\n",
      "    # Construct all type 3 matricies\n",
      "    for col in range(0, n):\n",
      "        for row in range(col + 1, n):\n",
      "            E = cmultadd(n, row, col, (-U[row, col] / U[col, col]))\n",
      "            E1 = cmultadd(n, row, col,  U[row, col] / U[col, col])\n",
      "            U = np.dot(E, U)\n",
      "            EL.append(E1)\n",
      "\n",
      "    # Construct all type 1 matrcies.\n",
      "    for j in range(0, n):\n",
      "        E = cmult(n, j, 1 / U[j, j])\n",
      "        E1 = cmult(n, j, U[j, j])\n",
      "        U = np.dot(E, U)\n",
      "        EL.append(E1)\n",
      "\n",
      "    for i in EL:\n",
      "        L = np.dot(L, i)\n",
      "\n",
      "    return [L, U]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A2 = np.random.randn(5, 5)\n",
      "L, U = LUDecomp(A2)\n",
      "np.allclose(A2, L.dot(U))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 46,
       "text": [
        "True"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Problem 3"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def detLU(mat):\n",
      "    \"\"\"\n",
      "    Use the LUDecomp function above to find the determinant of a matrix\n",
      "    \"\"\"\n",
      "    L, U = LUDecomp(mat)\n",
      "    sumL = 1\n",
      "    sumU = 1\n",
      "    for i in range(0, L.shape[0]):\n",
      "        sumL *= L[i, i]\n",
      "        sumU *= U[i ,i]\n",
      "    return sumL * sumU"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.allclose(detLU(A2), npla.det(A2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 48,
       "text": [
        "True"
       ]
      }
     ],
     "prompt_number": 48
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}